source: trunk/GSASIIstruct.py @ 379

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

finish spherical harmonics preferred orientation

  • Property svn:keywords set to Date Author Revision URL Id
File size: 77.5 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2011-09-22 15:11:18 +0000 (Thu, 22 Sep 2011) $
4# $Author: vondreele $
5# $Revision: 379 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 379 2011-09-22 15:11:18Z 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    def PrintSHPO(hapData):
559        print '\n Spherical harmonics preferred orientation: Order:' + \
560            str(hapData[4])+' Refine? '+str(hapData[2])
561        ptlbls = ' names :'
562        ptstr =  ' values:'
563        for item in hapData[5]:
564            ptlbls += '%12s'%(item)
565            ptstr += '%12.4f'%(hapData[5][item]) 
566        print ptlbls
567        print ptstr
568   
569    hapDict = {}
570    hapVary = []
571    controlDict = {}
572    poType = {}
573    poAxes = {}
574    spAxes = {}
575    spType = {}
576   
577    for phase in Phases:
578        HistoPhase = Phases[phase]['Histograms']
579        SGData = Phases[phase]['General']['SGData']
580        cell = Phases[phase]['General']['Cell'][1:7]
581        A = G2lat.cell2A(cell)
582        pId = Phases[phase]['pId']
583        for histogram in HistoPhase:
584            hapData = HistoPhase[histogram]
585            Histogram = Histograms[histogram]
586            hId = Histogram['hId']
587            limits = Histogram['Limits'][1]
588            inst = Histogram['Instrument Parameters']
589            inst = dict(zip(inst[3],inst[1]))
590            Zero = inst['Zero']
591            if 'C' in inst['Type']:
592                try:
593                    wave = inst['Lam']
594                except KeyError:
595                    wave = inst['Lam1']
596                dmin = wave/(2.0*sind(limits[1]/2.0))
597            pfx = str(pId)+':'+str(hId)+':'
598            for item in ['Scale','Extinction']:
599                hapDict[pfx+item] = hapData[item][0]
600                if hapData[item][1]:
601                    hapVary.append(pfx+item)
602            names = G2spc.HStrainNames(SGData)
603            for i,name in enumerate(names):
604                hapDict[pfx+name] = hapData['HStrain'][0][i]
605                if hapData['HStrain'][1][i]:
606                    hapVary.append(pfx+name)
607            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
608            if hapData['Pref.Ori.'][0] == 'MD':
609                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
610                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
611                if hapData['Pref.Ori.'][2]:
612                    hapVary.append(pfx+'MD')
613            else:                           #'SH' spherical harmonics
614                controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
615                controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
616                for item in hapData['Pref.Ori.'][5]:
617                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
618                    if hapData['Pref.Ori.'][2]:
619                        hapVary.append(pfx+item)
620            for item in ['Mustrain','Size']:
621                controlDict[pfx+item+'Type'] = hapData[item][0]
622                if hapData[item][0] in ['isotropic','uniaxial']:
623                    hapDict[pfx+item+':0'] = hapData[item][1][0]
624                    if hapData[item][2][0]:
625                        hapVary.append(pfx+item+':0')
626                    if hapData[item][0] == 'uniaxial':
627                        controlDict[pfx+item+'Axis'] = hapData[item][3]
628                        hapDict[pfx+item+':1'] = hapData[item][1][1]
629                        if hapData[item][2][1]:
630                            hapVary.append(pfx+item+':1')
631                else:       #generalized for mustrain or ellipsoidal for size
632                    if item == 'Mustrain':
633                        names = G2spc.MustrainNames(SGData)
634                        pwrs = []
635                        for name in names:
636                            h,k,l = name[1:]
637                            pwrs.append([int(h),int(k),int(l)])
638                        controlDict[pfx+'MuPwrs'] = pwrs
639                    for i in range(len(hapData[item][4])):
640                        sfx = ':'+str(i)
641                        hapDict[pfx+item+sfx] = hapData[item][4][i]
642                        if hapData[item][5][i]:
643                            hapVary.append(pfx+item+sfx)
644                           
645            if Print: 
646                print '\n Phase: ',phase,' in histogram: ',histogram
647                print '\n Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
648                print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
649                if hapData['Pref.Ori.'][0] == 'MD':
650                    Ax = hapData['Pref.Ori.'][3]
651                    print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
652                        ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
653                else: #'SH' for spherical harmonics
654                    PrintSHPO(hapData['Pref.Ori.'])
655                PrintSize(hapData['Size'])
656                PrintMuStrain(hapData['Mustrain'],SGData)
657                PrintHStrain(hapData['HStrain'],SGData)
658            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
659            refList = []
660            for h,k,l,d in HKLd:
661                ext,mul,Uniq = G2spc.GenHKLf([h,k,l],SGData)
662                if ext:
663                    continue
664                if 'C' in inst['Type']:
665                    pos = 2.0*asind(wave/(2.0*d))
666                    if limits[0] < pos < limits[1]:
667                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,Uniq])
668                else:
669                    raise ValueError 
670            Histogram['Reflection Lists'][phase] = refList
671    return hapVary,hapDict,controlDict
672   
673def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms):
674   
675    def PrintSizeAndSig(hapData,sizeSig):
676        line = '\n Size model:     '+hapData[0]
677        if hapData[0] in ['isotropic','uniaxial']:
678            line += ' equatorial:%12.3f'%(hapData[1][0])
679            if sizeSig[0][0]:
680                line += ', sig: %8.3f'%(sizeSig[0][0])
681            if hapData[0] == 'uniaxial':
682                line += ' axial:%12.3f'%(hapData[1][1])
683                if sizeSig[0][1]:
684                    line += ', sig: %8.3f'%(sizeSig[0][1])
685            print line
686        else:
687            print line
688            Snames = ['S11','S22','S33','S12','S13','S23']
689            ptlbls = ' name  :'
690            ptstr =  ' value :'
691            sigstr = ' sig   :'
692            for i,name in enumerate(Snames):
693                ptlbls += '%12s' % (name)
694                ptstr += '%12.6f' % (hapData[4][i])
695                if sizeSig[1][i]:
696                    sigstr += '%12.6f' % (sizeSig[1][i])
697                else:
698                    sigstr += 12*' '
699            print ptlbls
700            print ptstr
701            print sigstr
702       
703    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
704        line = '\n Mustrain model: '+hapData[0]
705        if hapData[0] in ['isotropic','uniaxial']:
706            line += ' equatorial:%12.1f'%(hapData[1][0])
707            if mustrainSig[0][0]:
708                line += ', sig: %8.1f'%(mustrainSig[0][0])
709            if hapData[0] == 'uniaxial':
710                line += ' axial:%12.1f'%(hapData[1][1])
711                if mustrainSig[0][1]:
712                     line += ', sig: %8.1f'%(mustrainSig[0][1])
713            print line
714        else:
715            print line
716            Snames = G2spc.MustrainNames(SGData)
717            ptlbls = ' name  :'
718            ptstr =  ' value :'
719            sigstr = ' sig   :'
720            for i,name in enumerate(Snames):
721                ptlbls += '%12s' % (name)
722                ptstr += '%12.6f' % (hapData[4][i])
723                if mustrainSig[1][i]:
724                    sigstr += '%12.6f' % (mustrainSig[1][i])
725                else:
726                    sigstr += 12*' '
727            print ptlbls
728            print ptstr
729            print sigstr
730           
731    def PrintHStrainAndSig(hapData,strainSig,SGData):
732        print '\n Hydrostatic strain: '
733        Hsnames = G2spc.HStrainNames(SGData)
734        ptlbls = ' name  :'
735        ptstr =  ' value :'
736        sigstr = ' sig   :'
737        for i,name in enumerate(Hsnames):
738            ptlbls += '%12s' % (name)
739            ptstr += '%12.6g' % (hapData[0][i])
740            if name in strainSig:
741                sigstr += '%12.6g' % (strainSig[name])
742            else:
743                sigstr += 12*' '
744        print ptlbls
745        print ptstr
746        print sigstr
747       
748    def PrintSHPOAndSig(hapData,POsig):
749        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
750        ptlbls = ' names :'
751        ptstr =  ' values:'
752        sigstr = ' sig   :'
753        for item in hapData[5]:
754            ptlbls += '%12s'%(item)
755            ptstr += '%12.4f'%(hapData[5][item])
756            if item in POsig:
757                sigstr += '%12.4f'%(POsig[item])
758            else:
759                sigstr += 12*' ' 
760        print ptlbls
761        print ptstr
762        print sigstr
763   
764    for phase in Phases:
765        HistoPhase = Phases[phase]['Histograms']
766        SGData = Phases[phase]['General']['SGData']
767        pId = Phases[phase]['pId']
768        for histogram in HistoPhase:
769            print '\n Phase: ',phase,' in histogram: ',histogram
770            hapData = HistoPhase[histogram]
771            Histogram = Histograms[histogram]
772            hId = Histogram['hId']
773            pfx = str(pId)+':'+str(hId)+':'
774           
775            PhFrExtPOSig = {}
776            for item in ['Scale','Extinction']:
777                hapData[item][0] = parmDict[pfx+item]
778                if pfx+item in sigDict:
779                    PhFrExtPOSig[item] = sigDict[pfx+item]
780            if hapData['Pref.Ori.'][0] == 'MD':
781                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
782                if pfx+'MD' in sigDict:
783                    PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
784            else:                           #'SH' spherical harmonics
785                for item in hapData['Pref.Ori.'][5]:
786                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
787                    if pfx+item in sigDict:
788                        PhFrExtPOSig[item] = sigDict[pfx+item]
789            print '\n'
790            if 'Scale' in PhFrExtPOSig:
791                print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
792            if 'Extinction' in PhFrExtPOSig:
793                print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
794            if hapData['Pref.Ori.'][0] == 'MD':
795                if 'MD' in PhFrExtPOSig:
796                    print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
797            else:
798                PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
799            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
800                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
801                'HStrain':{}}                 
802            for item in ['Mustrain','Size']:
803                if hapData[item][0] in ['isotropic','uniaxial']:                   
804                    hapData[item][1][0] = parmDict[pfx+item+':0']
805                    if item == 'Size':
806                        hapData[item][1][0] = min(10.,max(0.01,hapData[item][1][0]))
807                    if pfx+item+':0' in sigDict: 
808                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':0']
809                    if hapData[item][0] == 'uniaxial':
810                        hapData[item][1][1] = parmDict[pfx+item+':1']
811                        if item == 'Size':
812                            hapData[item][1][1] = min(10.,max(0.01,hapData[item][1][1]))                       
813                        if pfx+item+':1' in sigDict:
814                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':1']
815                else:       #generalized for mustrain or ellipsoidal for size
816                    for i in range(len(hapData[item][4])):
817                        sfx = ':'+str(i)
818                        hapData[item][4][i] = parmDict[pfx+item+sfx]
819                        if pfx+item+sfx in sigDict:
820                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
821            names = G2spc.HStrainNames(SGData)
822            for i,name in enumerate(names):
823                hapData['HStrain'][0][i] = parmDict[pfx+name]
824                if pfx+name in sigDict:
825                    SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
826            PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
827            PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
828            PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
829   
830def GetHistogramData(Histograms,Print=True):
831   
832    def GetBackgroundParms(hId,Background):
833        bakType,bakFlag = Background[:2]
834        backVals = Background[3:]
835        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
836        if bakFlag:                                 #returns backNames as varyList = backNames
837            return bakType,dict(zip(backNames,backVals)),backNames
838        else:                                       #no background varied; varyList = []
839            return bakType,dict(zip(backNames,backVals)),[]
840       
841    def GetInstParms(hId,Inst):
842        insVals,insFlags,insNames = Inst[1:4]
843        dataType = insVals[0]
844        instDict = {}
845        insVary = []
846        pfx = ':'+str(hId)+':'
847        for i,flag in enumerate(insFlags):
848            insName = pfx+insNames[i]
849            instDict[insName] = insVals[i]
850            if flag:
851                insVary.append(insName)
852        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.01)
853        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.01)
854        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
855        return dataType,instDict,insVary
856       
857    def GetSampleParms(hId,Sample):
858        sampVary = []
859        hfx = ':'+str(hId)+':'       
860        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius']}
861        Type = Sample['Type']
862        if 'Bragg' in Type:             #Bragg-Brentano
863            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
864                sampDict[hfx+item] = Sample[item][0]
865                if Sample[item][1]:
866                    sampVary.append(hfx+item)
867        elif 'Debye' in Type:        #Debye-Scherrer
868            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
869                sampDict[hfx+item] = Sample[item][0]
870                if Sample[item][1]:
871                    sampVary.append(hfx+item)
872        return Type,sampDict,sampVary
873       
874    def PrintBackground(Background):
875        print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
876        line = ' Coefficients: '
877        for i,back in enumerate(Background[3:]):
878            line += '%10.3f'%(back)
879            if i and not i%10:
880                line += '\n'+15*' '
881        print line
882       
883    def PrintInstParms(Inst):
884        print '\n Instrument Parameters:'
885        ptlbls = ' name  :'
886        ptstr =  ' value :'
887        varstr = ' refine:'
888        instNames = Inst[3][1:]
889        for i,name in enumerate(instNames):
890            ptlbls += '%12s' % (name)
891            ptstr += '%12.6f' % (Inst[1][i+1])
892            if name in ['Lam1','Lam2','Azimuth']:
893                varstr += 12*' '
894            else:
895                varstr += '%12s' % (str(bool(Inst[2][i+1])))
896        print ptlbls
897        print ptstr
898        print varstr
899       
900    def PrintSampleParms(Sample):
901        print '\n Sample Parameters:'
902        ptlbls = ' name  :'
903        ptstr =  ' value :'
904        varstr = ' refine:'
905        if 'Bragg' in Sample['Type']:
906            for item in ['Scale','Shift','Transparency']:
907                ptlbls += '%14s'%(item)
908                ptstr += '%14.4f'%(Sample[item][0])
909                varstr += '%14s'%(str(bool(Sample[item][1])))
910           
911        elif 'Debye' in Type:        #Debye-Scherrer
912            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
913                ptlbls += '%14s'%(item)
914                ptstr += '%14.4f'%(Sample[item][0])
915                varstr += '%14s'%(str(bool(Sample[item][1])))
916
917        print ptlbls
918        print ptstr
919        print varstr
920       
921
922    histDict = {}
923    histVary = []
924    controlDict = {}
925    for histogram in Histograms:
926        Histogram = Histograms[histogram]
927        hId = Histogram['hId']
928        pfx = ':'+str(hId)+':'
929        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
930       
931        Background = Histogram['Background'][0]
932        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
933        controlDict[pfx+'bakType'] = Type
934        histDict.update(bakDict)
935        histVary += bakVary
936       
937        Inst = Histogram['Instrument Parameters']
938        Type,instDict,insVary = GetInstParms(hId,Inst)
939        controlDict[pfx+'histType'] = Type
940        histDict.update(instDict)
941        histVary += insVary
942       
943        Sample = Histogram['Sample Parameters']
944        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
945        controlDict[pfx+'instType'] = Type
946        histDict.update(sampDict)
947        histVary += sampVary
948
949        if Print: 
950            print '\n Histogram: ',histogram,' histogram Id: ',hId
951            print 135*'-'
952            Units = {'C':' deg','T':' msec'}
953            units = Units[controlDict[pfx+'histType'][2]]
954            Limits = controlDict[pfx+'Limits']
955            print ' Instrument type: ',Sample['Type']
956            print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
957            PrintSampleParms(Sample)
958            PrintInstParms(Inst)
959            PrintBackground(Background)
960       
961    return histVary,histDict,controlDict
962   
963def SetHistogramData(parmDict,sigDict,Histograms):
964   
965    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
966        lenBack = len(Background[3:])
967        backSig = [0 for i in range(lenBack)]
968        for i in range(lenBack):
969            Background[3+i] = parmDict[pfx+'Back:'+str(i)]
970            if pfx+'Back:'+str(i) in sigDict:
971                backSig[i] = sigDict[pfx+'Back:'+str(i)]
972        return backSig
973       
974    def SetInstParms(pfx,Inst,parmDict,sigDict):
975        insVals,insFlags,insNames = Inst[1:4]
976        instSig = [0 for i in range(len(insVals))]
977        for i,flag in enumerate(insFlags):
978            insName = pfx+insNames[i]
979            insVals[i] = parmDict[insName]
980            if insName in sigDict:
981                instSig[i] = sigDict[insName]
982        return instSig
983       
984    def SetSampleParms(pfx,Sample,parmDict,sigDict):
985        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
986            sampSig = [0 for i in range(3)]
987            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
988                Sample[item][0] = parmDict[pfx+item]
989                if pfx+item in sigDict:
990                    sampSig[i] = sigDict[pfx+item]
991        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
992            sampSig = [0 for i in range(4)]
993            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
994                Sample[item][0] = parmDict[pfx+item]
995                if pfx+item in sigDict:
996                    sampSig[i] = sigDict[pfx+item]
997        return sampSig
998       
999    def PrintBackgroundSig(Background,backSig):
1000        print '\n Background function: ',Background[0]
1001        valstr = ' value : '
1002        sigstr = ' sig   : '
1003        for i,back in enumerate(Background[3:]):
1004            valstr += '%10.4f'%(back)
1005            if Background[1]:
1006                sigstr += '%10.4f'%(backSig[i])
1007            else:
1008                sigstr += 10*' '
1009        print valstr
1010        print sigstr
1011       
1012    def PrintInstParmsSig(Inst,instSig):
1013        print '\n Instrument Parameters:'
1014        ptlbls = ' names :'
1015        ptstr =  ' value :'
1016        sigstr = ' sig   :'
1017        instNames = Inst[3][1:]
1018        for i,name in enumerate(instNames):
1019            ptlbls += '%12s' % (name)
1020            ptstr += '%12.6f' % (Inst[1][i+1])
1021            if instSig[i+1]:
1022                sigstr += '%12.6f' % (instSig[i+1])
1023            else:
1024                sigstr += 12*' '
1025        print ptlbls
1026        print ptstr
1027        print sigstr
1028       
1029    def PrintSampleParmsSig(Sample,sampleSig):
1030        print '\n Sample Parameters:'
1031        ptlbls = ' names :'
1032        ptstr =  ' values:'
1033        sigstr = ' sig   :'
1034        if 'Bragg' in Sample['Type']:
1035            for i,item in enumerate(['Scale','Shift','Transparency']):
1036                ptlbls += '%14s'%(item)
1037                ptstr += '%14.4f'%(Sample[item][0])
1038                sigstr += '%14.4f'%(sampleSig[i])
1039           
1040        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1041            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1042                ptlbls += '%14s'%(item)
1043                ptstr += '%14.4f'%(Sample[item][0])
1044                sigstr += '%14.4f'%(sampleSig[i])
1045
1046        print ptlbls
1047        print ptstr
1048        print sigstr
1049       
1050    for histogram in Histograms:
1051        if 'PWDR' in histogram:
1052            Histogram = Histograms[histogram]
1053            hId = Histogram['hId']
1054            pfx = ':'+str(hId)+':'
1055            Background = Histogram['Background'][0]
1056            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1057           
1058            Inst = Histogram['Instrument Parameters']
1059            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1060       
1061            Sample = Histogram['Sample Parameters']
1062            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1063
1064            print '\n Histogram: ',histogram,' histogram Id: ',hId
1065            print 135*'-'
1066            print ' Instrument type: ',Sample['Type']
1067            PrintSampleParmsSig(Sample,sampSig)
1068            PrintInstParmsSig(Inst,instSig)
1069            PrintBackgroundSig(Background,backSig)
1070
1071
1072def SetupSFcalc(General,Atoms):
1073    ''' setup data for use in StructureFactor; mostly rearranging arrays
1074    input:
1075        General = dictionary of general phase info.
1076        Atoms = list of atom parameters
1077    returns:
1078        G = reciprocal metric tensor
1079        StrData = list of arrays; one entry per atom:
1080            T = atom types
1081            R = refinement flags, e.g. 'FXU'
1082            F = site fractions
1083            X = atom coordinates as numpy array
1084            M = atom multiplicities
1085            IA = isotropic/anisotropic thermal flags
1086            Uiso = isotropic thermal parameters if IA = 'I'; else = 0
1087            Uij = numpy array of 6 anisotropic thermal parameters if IA='A'; else = zeros
1088    '''           
1089    SGData = General['SGData']
1090    Cell = General['Cell']
1091    G,g = G2lat.cell2Gmat(Cell[1:7])        #skip refine & volume; get recip & real metric tensors
1092    cx,ct,cs,cia = General['AtomPtrs']
1093    X = [];F = [];T = [];IA = [];Uiso = [];Uij = [];R = [];M = []
1094    for atom in Atoms:
1095        T.append(atom[ct])
1096        R.append(atom[ct+1])
1097        F.append(atom[cx+3])
1098        X.append(np.array(atom[cx:cx+3]))
1099        M.append(atom[cia-1])
1100        IA.append(atom[cia])
1101        Uiso.append(atom[cia+1])
1102        Uij.append(np.array(atom[cia+2:cia+8]))
1103    return G,[T,R,np.array(F),np.array(X),np.array(M),IA,np.array(Uiso),np.array(Uij)]
1104           
1105def StructureFactor(H,G,SGData,StrData,FFtable):
1106    ''' Compute structure factor for a single h,k,l
1107    H = np.array(h,k,l)
1108    G = reciprocal metric tensor
1109    SGData = space group info. dictionary output from SpcGroup
1110    StrData = [
1111        [atom types],
1112        [refinement flags],
1113        [site fractions],
1114        np.array(coordinates),
1115        [multiplicities],
1116        [I/A flag],
1117        [Uiso],
1118        np.array(Uij)]
1119    FFtable = dictionary of form factor coeff. for atom types used in StrData
1120    '''
1121    twopi = 2.0*math.pi
1122    twopisq = 2.0*math.pi**2
1123    SQ = G2lat.calc_rDsq2(H,G)/4.0          # SQ = (sin(theta)/lambda)**2
1124    SQfactor = 8.0*SQ*math.pi**2
1125    print 'SQ',SQfactor
1126    FF = {}
1127    for type in FFtable.keys():
1128        FF[type] = G2el.ScatFac(FFtable[type],SQ)
1129    print 'FF',FF
1130    iabsnt,mulp,Uniq,Phs = G2spc.GenHKL(H,SGData,False)
1131    fa = [0,0,0]        #real
1132    fb = [0,0,0]        #imaginary
1133    if not iabsnt:
1134        phase = twopi*np.inner(Uniq,StrData[3])       #2pi[hx+ky+lz] for each atom in each equiv. hkl
1135        sinp = np.sin(phase)
1136        cosp = np.cos(phase)
1137        occ = StrData[2]*StrData[4]/mulp
1138        Tiso = np.multiply(StrData[6],SQfactor)
1139        Tuij = np.multiply(-SQfactor,np.inner(H,np.inner(G2spc.Uij2U(StrData[7]),H)))
1140        print 'sinp',sinp
1141        print 'cosp',cosp
1142        print 'occ',occ
1143        print 'Tiso',Tiso
1144        print 'Tuij',Tuij
1145    else:
1146        print 'Absent'
1147       
1148def Dict2Values(parmdict, varylist):
1149    '''Use before call to leastsq to setup list of values for the parameters
1150    in parmdict, as selected by key in varylist'''
1151    return [parmdict[key] for key in varylist] 
1152   
1153def Values2Dict(parmdict, varylist, values):
1154    ''' Use after call to leastsq to update the parameter dictionary with
1155    values corresponding to keys in varylist'''
1156    parmdict.update(zip(varylist,values))
1157   
1158def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1159    odfCor = 1.0
1160    H = refl[:3]
1161    cell = G2lat.Gmat2cell(g)
1162    Sangl = [0.,0.,0.]
1163    if 'Bragg' in calcControls[hfx+'instType']:
1164        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1165        IFCoup = True
1166    else:
1167        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1168        IFCoup = False
1169    phi,beta = G2lat.CrsAng(H,cell,SGData)
1170    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1171    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],None,calcControls[phfx+'SHord'])
1172    for item in SHnames:
1173        L,N = eval(item.strip('C'))
1174        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1175        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1176    return odfCor
1177   
1178def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1179    FORPI = 12.5663706143592
1180    odfCor = 1.0
1181    dFdODF = {}
1182    H = refl[:3]
1183    cell = G2lat.Gmat2cell(g)
1184    Sangl = [0.,0.,0.]
1185    if 'Bragg' in calcControls[hfx+'instType']:
1186        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1187        IFCoup = True
1188    else:
1189        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1190        IFCoup = False
1191    phi,beta = G2lat.CrsAng(H,cell,SGData)
1192    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1193    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],None,calcControls[phfx+'SHord'])
1194    for item in SHnames:
1195        L,N = eval(item.strip('C'))
1196        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1197        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1198        dFdODF[phfx+item] = FORPI*Kcsl
1199    return odfCor,dFdODF
1200   
1201def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1202    if calcControls[phfx+'poType'] == 'MD':
1203        MD = parmDict[phfx+'MD']
1204        MDAxis = calcControls[phfx+'MDAxis']
1205        sumMD = 0
1206        for H in refl[10]:           
1207            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1208            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1209            sumMD += A**3
1210        POcorr = sumMD/len(refl[10])
1211    else:   #spherical harmonics
1212        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1213    return POcorr
1214   
1215def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1216    POderv = {}
1217    if calcControls[phfx+'poType'] == 'MD':
1218        MD = parmDict[phfx+'MD']
1219        MDAxis = calcControls[phfx+'MDAxis']
1220        sumMD = 0
1221        sumdMD = 0
1222        for H in refl[10]:           
1223            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1224            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1225            sumMD += A**3
1226            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1227        POcorr = sumMD/len(refl[10])
1228        POderv[phfx+'MD'] = sumdMD/len(refl[10])
1229    else:   #spherical harmonics
1230        POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1231    return POcorr,POderv
1232   
1233def GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1234    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1235    Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1236    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1237    return Icorr
1238   
1239def GetIntensityDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1240    Icorr = GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1241    pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1242    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1243    dIdPola /= pola
1244    for iPO in dIdPO:
1245        dIdPO[iPO] /= POcorr
1246    dIdsh = 1./parmDict[hfx+'Scale']
1247    dIdsp = 1./parmDict[phfx+'Scale']
1248    return Icorr,dIdsh,dIdsp,dIdPola,dIdPO
1249       
1250def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1251    costh = cosd(refl[5]/2.)
1252    #crystallite size
1253    if calcControls[phfx+'SizeType'] == 'isotropic':
1254        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh)
1255    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1256        H = np.array(refl[:3])
1257        P = np.array(calcControls[phfx+'SizeAxis'])
1258        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1259        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']*costh)
1260        gam *= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)
1261    else:           #ellipsoidal crystallites - wrong not make sense
1262        H = np.array(refl[:3])
1263        gam += 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1264    #microstrain               
1265    if calcControls[phfx+'MustrainType'] == 'isotropic':
1266        gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
1267    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1268        H = np.array(refl[:3])
1269        P = np.array(calcControls[phfx+'MustrainAxis'])
1270        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1271        Si = parmDict[phfx+'Mustrain:0']
1272        Sa = parmDict[phfx+'Mustrain:1']
1273        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1274    else:       #generalized - P.W. Stephens model
1275        pwrs = calcControls[phfx+'MuPwrs']
1276        sum = 0
1277        for i,pwr in enumerate(pwrs):
1278            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1279        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
1280    return gam
1281       
1282def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1283    gamDict = {}
1284    costh = cosd(refl[5]/2.)
1285    tanth = tand(refl[5]/2.)
1286    #crystallite size derivatives
1287    if calcControls[phfx+'SizeType'] == 'isotropic':
1288        gamDict[phfx+'Size:0'] = -1.80*wave/(np.pi*costh)
1289    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1290        H = np.array(refl[:3])
1291        P = np.array(calcControls[phfx+'SizeAxis'])
1292        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1293        Si = parmDict[phfx+'Size:0']
1294        Sa = parmDict[phfx+'Size:1']
1295        gami = (1.80*wave/np.pi)/(Si*Sa)
1296        sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2)
1297        gam = gami*sqtrm/costh           
1298        gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)-gam/Si
1299        gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)-gam/Sa         
1300    else:           #ellipsoidal crystallites - do numerically? - not right not make sense
1301        H = np.array(refl[:3])
1302        gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1303    #microstrain derivatives               
1304    if calcControls[phfx+'MustrainType'] == 'isotropic':
1305        gamDict[phfx+'Mustrain:0'] =  0.018*tanth/np.pi           
1306    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1307        H = np.array(refl[:3])
1308        P = np.array(calcControls[phfx+'MustrainAxis'])
1309        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1310        Si = parmDict[phfx+'Mustrain:0']
1311        Sa = parmDict[phfx+'Mustrain:1']
1312        gami = 0.018*Si*Sa*tanth/np.pi
1313        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1314        gam = gami/sqtrm
1315        gamDict[phfx+'Mustrain:0'] = gam/Si-gami*Si*cosP**2/sqtrm**3
1316        gamDict[phfx+'Mustrain:1'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
1317    else:       #generalized - P.W. Stephens model
1318        pwrs = calcControls[phfx+'MuPwrs']
1319        const = 0.018*refl[4]**2*tanth
1320        for i,pwr in enumerate(pwrs):
1321            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1322    return gamDict
1323       
1324def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1325    h,k,l = refl[:3]
1326    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1327    d = np.sqrt(dsq)
1328    refl[4] = d
1329    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1330    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1331    if 'Bragg' in calcControls[hfx+'instType']:
1332        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1333            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1334    else:               #Debye-Scherrer - simple but maybe not right
1335        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1336    return pos
1337
1338def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1339    dpr = 180./np.pi
1340    h,k,l = refl[:3]
1341    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1342    dst = np.sqrt(dstsq)
1343    pos = refl[5]
1344    const = dpr/np.sqrt(1.0-wave*dst/4.0)
1345    dpdw = const*dst
1346    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1347    dpdA *= const*wave/(2.0*dst)
1348    dpdZ = 1.0
1349    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1350    if 'Bragg' in calcControls[hfx+'instType']:
1351        dpdSh = -4.*const*cosd(pos/2.0)
1352        dpdTr = -const*sind(pos)*100.0
1353        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1354    else:               #Debye-Scherrer - simple but maybe not right
1355        dpdXd = -const*cosd(pos)
1356        dpdYd = -const*sind(pos)
1357        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1358           
1359def GetHStrainShift(refl,SGData,phfx,parmDict):
1360    laue = SGData['SGLaue']
1361    uniq = SGData['SGUniq']
1362    h,k,l = refl[:3]
1363    if laue in ['m3','m3m']:
1364        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)
1365    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1366        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1367    elif laue in ['3R','3mR']:
1368        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1369    elif laue in ['4/m','4/mmm']:
1370        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1371    elif laue in ['mmm']:
1372        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1373    elif laue in ['2/m']:
1374        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1375        if uniq == 'a':
1376            Dij += parmDict[phfx+'D23']*k*l
1377        elif uniq == 'b':
1378            Dij += parmDict[phfx+'D13']*h*l
1379        elif uniq == 'c':
1380            Dij += parmDict[phfx+'D12']*h*k
1381    else:
1382        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1383            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1384    return Dij*refl[4]**2*tand(refl[5]/2.0)
1385           
1386def GetHStrainShiftDerv(refl,SGData,phfx):
1387    laue = SGData['SGLaue']
1388    uniq = SGData['SGUniq']
1389    h,k,l = refl[:3]
1390    if laue in ['m3','m3m']:
1391        dDijDict = {phfx+'D11':h**2+k**2+l**2,}
1392    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1393        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1394    elif laue in ['3R','3mR']:
1395        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1396    elif laue in ['4/m','4/mmm']:
1397        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1398    elif laue in ['mmm']:
1399        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1400    elif laue in ['2/m']:
1401        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1402        if uniq == 'a':
1403            dDijDict[phfx+'D23'] = k*l
1404        elif uniq == 'b':
1405            dDijDict[phfx+'D13'] = h*l
1406        elif uniq == 'c':
1407            dDijDict[phfx+'D12'] = h*k
1408            names.append()
1409    else:
1410        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1411            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1412    for item in dDijDict:
1413        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
1414    return dDijDict
1415           
1416def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1417   
1418    def GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse):
1419        U = parmDict[hfx+'U']
1420        V = parmDict[hfx+'V']
1421        W = parmDict[hfx+'W']
1422        X = parmDict[hfx+'X']
1423        Y = parmDict[hfx+'Y']
1424        tanPos = tand(refl[5]/2.0)
1425        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
1426        sig = max(0.001,sig)
1427        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
1428        gam = max(0.01,gam)
1429        return sig,gam
1430               
1431    hId = Histogram['hId']
1432    hfx = ':%d:'%(hId)
1433    bakType = calcControls[hfx+'bakType']
1434    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1435    yc = np.zeros_like(yb)
1436       
1437    if 'C' in calcControls[hfx+'histType']:   
1438        shl = max(parmDict[hfx+'SH/L'],0.0005)
1439        Ka2 = False
1440        if hfx+'Lam1' in parmDict.keys():
1441            wave = parmDict[hfx+'Lam1']
1442            Ka2 = True
1443            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1444            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1445        else:
1446            wave = parmDict[hfx+'Lam']
1447    else:
1448        print 'TOF Undefined at present'
1449        raise ValueError
1450    for phase in Histogram['Reflection Lists']:
1451        refList = Histogram['Reflection Lists'][phase]
1452        Phase = Phases[phase]
1453        pId = Phase['pId']
1454        pfx = '%d::'%(pId)
1455        phfx = '%d:%d:'%(pId,hId)
1456        SGData = Phase['General']['SGData']
1457        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1458        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1459        sizeEllipse = []
1460        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1461            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1462        for refl in refList:
1463            if 'C' in calcControls[hfx+'histType']:
1464                h,k,l = refl[:3]
1465                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1466                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1467                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
1468                Icorr = GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1469                if 'Pawley' in Phase['General']['Type']:
1470                    try:
1471                        refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
1472                    except KeyError:
1473#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1474                        continue
1475                else:
1476                    raise ValueError       #wants strctrfacr calc here
1477                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1478                iBeg = np.searchsorted(x,refl[5]-fmin)
1479                iFin = np.searchsorted(x,refl[5]+fmax)
1480                if not iBeg+iFin:       #peak below low limit - skip peak
1481                    continue
1482                elif not iBeg-iFin:     #peak above high limit - done
1483                    break
1484                yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1485                if Ka2:
1486                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1487                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
1488                    iBeg = np.searchsorted(x,pos2-fmin)
1489                    iFin = np.searchsorted(x,pos2+fmax)
1490                    if not iBeg+iFin:       #peak below low limit - skip peak
1491                        continue
1492                    elif not iBeg-iFin:     #peak above high limit - done
1493                        return yc,yb
1494                    yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
1495            elif 'T' in calcControls[hfx+'histType']:
1496                print 'TOF Undefined at present'
1497                raise ValueError    #no TOF yet
1498    return yc,yb   
1499           
1500def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1501   
1502    def cellVaryDerv(pfx,SGData,dpdA): 
1503        if SGData['SGLaue'] in ['-1',]:
1504            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1505                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1506        elif SGData['SGLaue'] in ['2/m',]:
1507            if SGData['SGUniq'] == 'a':
1508                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1509            elif SGData['SGUniq'] == 'b':
1510                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1511            else:
1512                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1513        elif SGData['SGLaue'] in ['mmm',]:
1514            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1515        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1516            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
1517        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1518            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
1519        elif SGData['SGLaue'] in ['3R', '3mR']:
1520            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1521        elif SGData['SGLaue'] in ['m3m','m3']:
1522            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
1523   
1524    lenX = len(x)               
1525    hId = Histogram['hId']
1526    hfx = ':%d:'%(hId)
1527    bakType = calcControls[hfx+'bakType']
1528    dMdv = np.zeros(shape=(len(varylist),len(x)))
1529    if hfx+'Back:0' in varylist:
1530        dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1531        bBpos =varylist.index(hfx+'Back:0')
1532        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1533       
1534    if 'C' in calcControls[hfx+'histType']:   
1535        dx = x[1]-x[0]
1536        shl = max(parmDict[hfx+'SH/L'],0.002)
1537        Ka2 = False
1538        if hfx+'Lam1' in parmDict.keys():
1539            wave = parmDict[hfx+'Lam1']
1540            Ka2 = True
1541            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1542            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1543        else:
1544            wave = parmDict[hfx+'Lam']
1545    else:
1546        print 'TOF Undefined at present'
1547        raise ValueError
1548    for phase in Histogram['Reflection Lists']:
1549        refList = Histogram['Reflection Lists'][phase]
1550        Phase = Phases[phase]
1551        SGData = Phase['General']['SGData']
1552        pId = Phase['pId']
1553        pfx = '%d::'%(pId)
1554        phfx = '%d:%d:'%(pId,hId)
1555        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1556        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1557        sizeEllipse = []
1558        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1559            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1560        for refl in refList:
1561            if 'C' in calcControls[hfx+'histType']:        #CW powder
1562                h,k,l = refl[:3]
1563                iref = pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]
1564                Icorr,dIdsh,dIdsp,dIdpola,dIdPO = GetIntensityDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1565                if 'Pawley' in Phase['General']['Type']:
1566                    try:
1567                        refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
1568                    except KeyError:
1569#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1570                        continue
1571                else:
1572                    raise ValueError       #wants strctrfacr deriv calc here
1573                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1574                iBeg = np.searchsorted(x,refl[5]-fmin)
1575                iFin = np.searchsorted(x,refl[5]+fmax)
1576                if not iBeg+iFin:       #peak below low limit - skip peak
1577                    continue
1578                elif not iBeg-iFin:     #peak above high limit - done
1579                    break
1580                pos = refl[5]
1581                tanth = tand(pos/2.0)
1582                costh = cosd(pos/2.0)
1583                dMdpk = np.zeros(shape=(6,len(x)))
1584                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
1585                for i in range(1,5):
1586                    dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i]
1587                dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[0]
1588                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1589                if Ka2:
1590                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1591                    kdelt = int((pos2-refl[5])/dx)               
1592                    iBeg2 = min(lenX,iBeg+kdelt)
1593                    iFin2 = min(lenX,iFin+kdelt)
1594                    if iBeg2-iFin2:
1595                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
1596                        for i in range(1,5):
1597                            dMdpk[i][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]
1598                        dMdpk[0][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[0]
1599                        dMdpk[5][iBeg2:iFin2] += 100.*dx*Icorr*dMdipk2[0]
1600                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[8]}
1601                try:
1602                    idx = varylist.index(pfx+'PWLref:'+str(iref))
1603                    dMdv[idx] = dervDict['int']/refl[8]
1604                except ValueError:
1605                    pass
1606                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1607                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1608                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1609                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1610                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1611                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1612                    hfx+'DisplaceY':[dpdY,'pos'],}
1613                for name in names:
1614                    if name in varylist:
1615                        item = names[name]
1616                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
1617                for iPO in dIdPO:
1618                    if iPO in varylist:
1619                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
1620                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1621                for name,dpdA in cellDervNames:
1622                    if name in varylist:
1623                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
1624                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1625                for name in dDijDict:
1626                    if name in varylist:
1627                        dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos']
1628                gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse)
1629                for name in gamDict:
1630                    if name in varylist:
1631                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']                       
1632            elif 'T' in calcControls[hfx+'histType']:
1633                print 'TOF Undefined at present'
1634                raise ValueError    #no TOF yet
1635           
1636    return dMdv   
1637                   
1638def Refine(GPXfile,dlg):
1639    import cPickle
1640    import pytexture as ptx
1641    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
1642   
1643    def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
1644        parmdict.update(zip(varylist,values))
1645        G2mv.Dict2Map(parmDict)
1646        Histograms,Phases = HistoPhases
1647        dMdv = np.empty(0)
1648        for histogram in Histograms:
1649            if 'PWDR' in histogram[:4]:
1650                Histogram = Histograms[histogram]
1651                hId = Histogram['hId']
1652                hfx = ':%d:'%(hId)
1653                Limits = calcControls[hfx+'Limits']
1654                x,y,w,yc,yb,yd = Histogram['Data']
1655                xB = np.searchsorted(x,Limits[0])
1656                xF = np.searchsorted(x,Limits[1])
1657                dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
1658                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1659                if len(dMdv):
1660                    dMdv = np.concatenate((dMdv,dMdvh))
1661                else:
1662                    dMdv = dMdvh
1663        return dMdv
1664   
1665    def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
1666        parmdict.update(zip(varylist,values))
1667        G2mv.Dict2Map(parmDict)
1668        Histograms,Phases = HistoPhases
1669        M = np.empty(0)
1670        sumwYo = 0
1671        Nobs = 0
1672        for histogram in Histograms:
1673            if 'PWDR' in histogram[:4]:
1674                Histogram = Histograms[histogram]
1675                hId = Histogram['hId']
1676                hfx = ':%d:'%(hId)
1677                Limits = calcControls[hfx+'Limits']
1678                x,y,w,yc,yb,yd = Histogram['Data']
1679                yc *= 0.0                           #zero full calcd profiles
1680                yb *= 0.0
1681                yd *= 0.0
1682                xB = np.searchsorted(x,Limits[0])
1683                xF = np.searchsorted(x,Limits[1])
1684                Nobs += xF-xB
1685                Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
1686                sumwYo += Histogram['sumwYo']
1687                yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
1688                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1689                yc[xB:xF] *= parmDict[hfx+'Scale']
1690                yc[xB:xF] += yb[xB:xF]
1691                yd[xB:xF] = yc[xB:xF]-y[xB:xF]          #yc-yo then all dydv have no '-' needed
1692                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
1693                M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF])))
1694        Histograms['sumwYo'] = sumwYo
1695        Histograms['Nobs'] = Nobs
1696        Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
1697        if dlg:
1698            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
1699            if not GoOn:
1700                return -M           #abort!!
1701        return M
1702   
1703    ShowBanner()
1704    varyList = []
1705    parmDict = {}
1706    calcControls = {}   
1707    Controls = GetControls(GPXfile)
1708    ShowControls(Controls)           
1709    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
1710    if not Phases:
1711        print ' *** ERROR - you have no histograms to refine! ***'
1712        print ' *** Refine aborted ***'
1713        raise Exception
1714    if not Histograms:
1715        print ' *** ERROR - you have no data to refine with! ***'
1716        print ' *** Refine aborted ***'
1717        raise Exception
1718    phaseVary,phaseDict,pawleyLookup = GetPhaseData(Phases)
1719    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
1720    calcControls.update(controlDict)
1721    histVary,histDict,controlDict = GetHistogramData(Histograms)
1722    calcControls.update(controlDict)
1723    varyList = phaseVary+hapVary+histVary
1724    parmDict.update(phaseDict)
1725    parmDict.update(hapDict)
1726    parmDict.update(histDict)
1727    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
1728    groups,parmlist = G2mv.GroupConstraints(constrDict)
1729    G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
1730    G2mv.Map2Dict(parmDict,varyList)
1731
1732    while True:
1733        begin = time.time()
1734        values =  np.array(Dict2Values(parmDict, varyList))
1735        Ftol = Controls['min dM/M']
1736        Factor = Controls['shift factor']
1737        if Controls['deriv type'] == 'analytic':
1738            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
1739                ftol=Ftol,col_deriv=True,factor=Factor,
1740                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1741            ncyc = int(result[2]['nfev']/2)               
1742        else:           #'numeric'
1743            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
1744                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1745            ncyc = int(result[2]['nfev']/len(varyList))
1746#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
1747#        for item in table: print item,table[item]               #useful debug - are things shifting?
1748        runtime = time.time()-begin
1749        chisq = np.sum(result[2]['fvec']**2)
1750        Values2Dict(parmDict, varyList, result[0])
1751        G2mv.Dict2Map(parmDict)
1752        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
1753        GOF = chisq/(Histograms['Nobs']-len(varyList))
1754        print '\n Refinement results:'
1755        print 135*'-'
1756        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
1757        print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1758        print ' Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1759        try:
1760            covMatrix = result[1]*GOF
1761            sig = np.sqrt(np.diag(covMatrix))
1762            xvar = np.outer(sig,np.ones_like(sig))
1763            cov = np.divide(np.divide(covMatrix,xvar),xvar.T)
1764            if np.any(np.isnan(sig)):
1765                print '*** Least squares aborted - some invalid esds possible ***'
1766#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
1767#            for item in table: print item,table[item]               #useful debug - are things shifting?
1768            break                   #refinement succeeded - finish up!
1769        except TypeError:          #result[1] is None on singular matrix
1770            print '**** Refinement failed - singular matrix ****'
1771            Ipvt = result[2]['ipvt']
1772            for i,ipvt in enumerate(Ipvt):
1773                if not np.sum(result[2]['fjac'],axis=1)[i]:
1774                    print 'Removing parameter: ',varyList[ipvt-1]
1775                    del(varyList[ipvt-1])
1776                    break
1777
1778#    print 'dependentParmList: ',G2mv.dependentParmList
1779#    print 'arrayList: ',G2mv.arrayList
1780#    print 'invarrayList: ',G2mv.invarrayList
1781#    print 'indParmList: ',G2mv.indParmList
1782#    print 'fixedDict: ',G2mv.fixedDict
1783    sigDict = dict(zip(varyList,sig))
1784    SetPhaseData(parmDict,sigDict,Phases)
1785    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
1786    SetHistogramData(parmDict,sigDict,Histograms)
1787    covData = {'variables':result[0],'varyList':varyList,'covariance':cov}
1788    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
1789#for testing purposes!!!
1790#    file = open('structTestdata.dat','wb')
1791#    cPickle.dump(parmDict,file,1)
1792#    cPickle.dump(varyList,file,1)
1793#    for histogram in Histograms:
1794#        if 'PWDR' in histogram[:4]:
1795#            Histogram = Histograms[histogram]
1796#    cPickle.dump(Histogram,file,1)
1797#    cPickle.dump(Phases,file,1)
1798#    cPickle.dump(calcControls,file,1)
1799#    cPickle.dump(pawleyLookup,file,1)
1800#    file.close()
1801
1802def main():
1803    arg = sys.argv
1804    if len(arg) > 1:
1805        GPXfile = arg[1]
1806        if not ospath.exists(GPXfile):
1807            print 'ERROR - ',GPXfile," doesn't exist!"
1808            exit()
1809        GPXpath = ospath.dirname(arg[1])
1810        Refine(GPXfile)
1811    else:
1812        print 'ERROR - missing filename'
1813        exit()
1814    print "Done"
1815         
1816if __name__ == '__main__':
1817    main()
Note: See TracBrowser for help on using the repository browser.