source: trunk/GSASIIstruct.py @ 424

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

remove 'Back'ground from consideration for constraints
add a bit more for cell parms in seq refinements

  • Property svn:keywords set to Date Author Revision URL Id
File size: 120.1 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2011-11-18 20:31:18 +0000 (Fri, 18 Nov 2011) $
4# $Author: vondreele $
5# $Revision: 424 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 424 2011-11-18 20:31:18Z vondreele $
8########### SVN repository information ###################
9import sys
10import os
11import os.path as ospath
12import numpy as np
13import numpy.linalg as nl
14import cPickle
15import time
16import math
17import GSASIIpath
18import GSASIIElem as G2el
19import GSASIIlattice as G2lat
20import GSASIIspc as G2spc
21import GSASIIpwd as G2pwd
22import GSASIImapvars as G2mv
23import scipy.optimize as so
24
25sind = lambda x: np.sin(x*np.pi/180.)
26cosd = lambda x: np.cos(x*np.pi/180.)
27tand = lambda x: np.tan(x*np.pi/180.)
28asind = lambda x: 180.*np.arcsin(x)/np.pi
29atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
30
31
32def ShowBanner():
33    print 80*'*'
34    print '   General Structure Analysis System-II Crystal Structure Refinement'
35    print '     by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
36    print ' This product includes software developed by the UChicago Argonne, LLC,' 
37    print '            as Operator of Argonne National Laboratory.'
38    print 80*'*','\n'
39
40def GetControls(GPXfile):
41    ''' Returns dictionary of control items found in GSASII gpx file
42    input:
43        GPXfile = .gpx full file name
44    return:
45        Controls = dictionary of control items
46    '''
47    Controls = {'deriv type':'analytical','min dM/M':0.0001,'shift factor':1.}
48    file = open(GPXfile,'rb')
49    while True:
50        try:
51            data = cPickle.load(file)
52        except EOFError:
53            break
54        datum = data[0]
55        if datum[0] == 'Controls':
56            Controls.update(datum[1])
57    file.close()
58    return Controls
59   
60def ShowControls(Controls):
61    print ' Least squares controls:'
62    print ' Derivative type: ',Controls['deriv type']
63    print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
64    print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
65   
66def GetConstraints(GPXfile):
67    constList = []
68    file = open(GPXfile,'rb')
69    while True:
70        try:
71            data = cPickle.load(file)
72        except EOFError:
73            break
74        datum = data[0]
75        if datum[0] == 'Constraints':
76            constDict = datum[1]
77            for item in constDict:
78                constList += constDict[item]
79    file.close()
80    constDict = []
81    constFlag = []
82    fixedList = []
83    for item in constList:
84        if item[-2]:
85            fixedList.append(str(item[-2]))
86        else:
87            fixedList.append('0')
88        if item[-1]:
89            constFlag.append(['VARY'])
90        else:
91            constFlag.append([])
92        itemDict = {}
93        for term in item[:-2]:
94            itemDict[term[1]] = term[0]
95        constDict.append(itemDict)
96    return constDict,constFlag,fixedList
97   
98def GetPhaseNames(GPXfile):
99    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
100    input:
101        GPXfile = gpx full file name
102    return:
103        PhaseNames = list of phase names
104    '''
105    file = open(GPXfile,'rb')
106    PhaseNames = []
107    while True:
108        try:
109            data = cPickle.load(file)
110        except EOFError:
111            break
112        datum = data[0]
113        if 'Phases' == datum[0]:
114            for datus in data[1:]:
115                PhaseNames.append(datus[0])
116    file.close()
117    return PhaseNames
118
119def GetAllPhaseData(GPXfile,PhaseName):
120    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
121    input:
122        GPXfile = gpx full file name
123        PhaseName = phase name
124    return:
125        phase dictionary
126    '''       
127    file = open(GPXfile,'rb')
128    General = {}
129    Atoms = []
130    while True:
131        try:
132            data = cPickle.load(file)
133        except EOFError:
134            break
135        datum = data[0]
136        if 'Phases' == datum[0]:
137            for datus in data[1:]:
138                if datus[0] == PhaseName:
139                    break
140    file.close()
141    return datus[1]
142   
143def GetHistograms(GPXfile,hNames):
144    """ Returns a dictionary of histograms found in GSASII gpx file
145    input:
146        GPXfile = .gpx full file name
147        hNames = list of histogram names
148    return:
149        Histograms = dictionary of histograms (types = PWDR & HKLF)
150    """
151    file = open(GPXfile,'rb')
152    Histograms = {}
153    while True:
154        try:
155            data = cPickle.load(file)
156        except EOFError:
157            break
158        datum = data[0]
159        hist = datum[0]
160        if hist in hNames:
161            if 'PWDR' in hist[:4]:
162                PWDRdata = {}
163                PWDRdata['Data'] = datum[1][1]          #powder data arrays
164                PWDRdata[data[2][0]] = data[2][1]       #Limits
165                PWDRdata[data[3][0]] = data[3][1]       #Background
166                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
167                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
168                try:
169                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
170                except IndexError:
171                    PWDRdata['Reflection lists'] = {}
172   
173                Histograms[hist] = PWDRdata
174            elif 'HKLF' in hist[:4]:
175                HKLFdata = []
176                datum = data[0]
177                HKLFdata = datum[1:][0]
178                Histograms[hist] = HKLFdata           
179    file.close()
180    return Histograms
181   
182def GetHistogramNames(GPXfile,hType):
183    """ Returns a list of histogram names found in GSASII gpx file
184    input:
185        GPXfile = .gpx full file name
186        hType = list ['PWDR','HKLF']
187    return:
188        HistogramNames = list of histogram names (types = PWDR & HKLF)
189    """
190    file = open(GPXfile,'rb')
191    HistogramNames = []
192    while True:
193        try:
194            data = cPickle.load(file)
195        except EOFError:
196            break
197        datum = data[0]
198        if datum[0][:4] in hType:
199            HistogramNames.append(datum[0])
200    file.close()
201    return HistogramNames
202   
203def GetUsedHistogramsAndPhases(GPXfile):
204    ''' Returns all histograms that are found in any phase
205    and any phase that uses a histogram
206    input:
207        GPXfile = .gpx full file name
208    return:
209        Histograms = dictionary of histograms as {name:data,...}
210        Phases = dictionary of phases that use histograms
211    '''
212    phaseNames = GetPhaseNames(GPXfile)
213    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
214    allHistograms = GetHistograms(GPXfile,histoList)
215    phaseData = {}
216    for name in phaseNames: 
217        phaseData[name] =  GetAllPhaseData(GPXfile,name)
218    Histograms = {}
219    Phases = {}
220    for phase in phaseData:
221        Phase = phaseData[phase]
222        if Phase['Histograms']:
223            if phase not in Phases:
224                pId = phaseNames.index(phase)
225                Phase['pId'] = pId
226                Phases[phase] = Phase
227            for hist in Phase['Histograms']:
228                if hist not in Histograms:
229                    Histograms[hist] = allHistograms[hist]
230                    #future restraint, etc. histograms here           
231                    hId = histoList.index(hist)
232                    Histograms[hist]['hId'] = hId
233    return Histograms,Phases
234   
235def GPXBackup(GPXfile,makeBack=True):
236    import distutils.file_util as dfu
237    GPXpath,GPXname = ospath.split(GPXfile)
238    if GPXpath == '': GPXpath = '.'
239    Name = ospath.splitext(GPXname)[0]
240    files = os.listdir(GPXpath)
241    last = 0
242    for name in files:
243        name = name.split('.')
244        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
245            if makeBack:
246                last = max(last,int(name[1].strip('bak'))+1)
247            else:
248                last = max(last,int(name[1].strip('bak')))
249    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
250    dfu.copy_file(GPXfile,GPXback)
251    return GPXback
252       
253def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
254    ''' Updates gpxfile from all histograms that are found in any phase
255    and any phase that used a histogram
256    input:
257        GPXfile = .gpx full file name
258        Histograms = dictionary of histograms as {name:data,...}
259        Phases = dictionary of phases that use histograms
260        CovData = dictionary of refined variables, varyList, & covariance matrix
261        makeBack = True if new backup of .gpx file is to be made; else use the last one made
262    '''
263                       
264    GPXback = GPXBackup(GPXfile,makeBack)
265    print '\n',135*'-'
266    print 'Read from file:',GPXback
267    print 'Save to file  :',GPXfile
268    infile = open(GPXback,'rb')
269    outfile = open(GPXfile,'wb')
270    while True:
271        try:
272            data = cPickle.load(infile)
273        except EOFError:
274            break
275        datum = data[0]
276#        print 'read: ',datum[0]
277        if datum[0] == 'Phases':
278            for iphase in range(len(data)):
279                if data[iphase][0] in Phases:
280                    phaseName = data[iphase][0]
281                    data[iphase][1].update(Phases[phaseName])
282        elif datum[0] == 'Covariance':
283            data[0][1] = CovData
284        try:
285            histogram = Histograms[datum[0]]
286#            print 'found ',datum[0]
287            data[0][1][1] = histogram['Data']
288            for datus in data[1:]:
289#                print '    read: ',datus[0]
290                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
291                    datus[1] = histogram[datus[0]]
292        except KeyError:
293            pass
294                               
295        cPickle.dump(data,outfile,1)
296    infile.close()
297    outfile.close()
298    print 'refinement save successful'
299   
300def SetSeqResult(GPXfile,Histograms,SeqResult):
301    GPXback = GPXBackup(GPXfile)
302    print '\n',135*'-'
303    print 'Read from file:',GPXback
304    print 'Save to file  :',GPXfile
305    infile = open(GPXback,'rb')
306    outfile = open(GPXfile,'wb')
307    while True:
308        try:
309            data = cPickle.load(infile)
310        except EOFError:
311            break
312        datum = data[0]
313        if datum[0] == 'Sequental results':
314            data[0][1] = SeqResult
315        try:
316            histogram = Histograms[datum[0]]
317            data[0][1][1] = histogram['Data']
318            for datus in data[1:]:
319                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
320                    datus[1] = histogram[datus[0]]
321        except KeyError:
322            pass
323                               
324        cPickle.dump(data,outfile,1)
325    infile.close()
326    outfile.close()
327    print 'refinement save successful'
328   
329                   
330def GetPWDRdata(GPXfile,PWDRname):
331    ''' Returns powder data from GSASII gpx file
332    input:
333        GPXfile = .gpx full file name
334        PWDRname = powder histogram name as obtained from GetHistogramNames
335    return:
336        PWDRdata = powder data dictionary with:
337            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
338       
339    '''
340    file = open(GPXfile,'rb')
341    PWDRdata = {}
342    while True:
343        try:
344            data = cPickle.load(file)
345        except EOFError:
346            break
347        datum = data[0]
348        if datum[0] == PWDRname:
349            PWDRdata['Data'] = datum[1][1]          #powder data arrays
350            PWDRdata[data[2][0]] = data[2][1]       #Limits
351            PWDRdata[data[3][0]] = data[3][1]       #Background
352            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
353            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
354            try:
355                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
356            except IndexError:
357                PWDRdata['Reflection lists'] = {}
358    file.close()
359    return PWDRdata
360   
361def GetHKLFdata(GPXfile,HKLFname):
362    ''' Returns single crystal data from GSASII gpx file
363    input:
364        GPXfile = .gpx full file name
365        HKLFname = single crystal histogram name as obtained from GetHistogramNames
366    return:
367        HKLFdata = single crystal data list of reflections: for each reflection:
368            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
369    '''
370    file = open(GPXfile,'rb')
371    HKLFdata = []
372    while True:
373        try:
374            data = cPickle.load(file)
375        except EOFError:
376            break
377        datum = data[0]
378        if datum[0] == HKLFname:
379            HKLFdata = datum[1:][0]
380    file.close()
381    return HKLFdata
382   
383def GetFFtable(General):
384    ''' returns a dictionary of form factor data for atom types found in General
385    input:
386        General = dictionary of phase info.; includes AtomTypes
387    return:
388        FFtable = dictionary of form factor data; key is atom type
389    '''
390    atomTypes = General['AtomTypes']
391    FFtable = {}
392    for El in atomTypes:
393        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
394        for item in FFs:
395            if item['Symbol'] == El.upper():
396                FFtable[El] = item
397    return FFtable
398   
399def GetBLtable(General):
400    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
401    input:
402        General = dictionary of phase info.; includes AtomTypes & Isotopes
403    return:
404        BLtable = dictionary of scattering length data; key is atom type
405    '''
406    atomTypes = General['AtomTypes']
407    BLtable = {}
408    isotopes = General['Isotopes']
409    isotope = General['Isotope']
410    for El in atomTypes:
411        BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
412    return BLtable
413       
414def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
415    if SGLaue in ['-1','2/m','mmm']:
416        return                      #no Pawley symmetry required constraints
417    for i,varyI in enumerate(pawleyVary):
418        refI = int(varyI.split(':')[-1])
419        ih,ik,il = PawleyRef[refI][:3]
420        for varyJ in pawleyVary[0:i]:
421            refJ = int(varyJ.split(':')[-1])
422            jh,jk,jl = PawleyRef[refJ][:3]
423            if SGLaue in ['4/m','4/mmm']:
424                isum = ih**2+ik**2
425                jsum = jh**2+jk**2
426                if abs(il) == abs(jl) and isum == jsum:
427                    G2mv.StoreEquivalence(varyJ,(varyI,))
428            elif SGLaue in ['3R','3mR']:
429                isum = ih**2+ik**2+il**2
430                jsum = jh**2+jk**2*jl**2
431                isum2 = ih*ik+ih*il+ik*il
432                jsum2 = jh*jk+jh*jl+jk*jl
433                if isum == jsum and isum2 == jsum2:
434                    G2mv.StoreEquivalence(varyJ,(varyI,))
435            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
436                isum = ih**2+ik**2+ih*ik
437                jsum = jh**2+jk**2+jh*jk
438                if abs(il) == abs(jl) and isum == jsum:
439                    G2mv.StoreEquivalence(varyJ,(varyI,))
440            elif SGLaue in ['m3','m3m']:
441                isum = ih**2+ik**2+il**2
442                jsum = jh**2+jk**2+jl**2
443                if isum == jsum:
444                    G2mv.StoreEquivalence(varyJ,(varyI,))
445                   
446def cellVary(pfx,SGData): 
447    if SGData['SGLaue'] in ['-1',]:
448        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
449    elif SGData['SGLaue'] in ['2/m',]:
450        if SGData['SGUniq'] == 'a':
451            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
452        elif SGData['SGUniq'] == 'b':
453            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
454        else:
455            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
456    elif SGData['SGLaue'] in ['mmm',]:
457        return [pfx+'A0',pfx+'A1',pfx+'A2']
458    elif SGData['SGLaue'] in ['4/m','4/mmm']:
459        return [pfx+'A0',pfx+'A2']
460    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
461        return [pfx+'A0',pfx+'A2']
462    elif SGData['SGLaue'] in ['3R', '3mR']:
463        return [pfx+'A0',pfx+'A3']                       
464    elif SGData['SGLaue'] in ['m3m','m3']:
465        return [pfx+'A0',]
466       
467           
468def GetPhaseData(PhaseData,Print=True):
469           
470    def PrintFFtable(FFtable):
471        print '\n X-ray scattering factors:'
472        print '   Symbol     fa                                      fb                                      fc'
473        print 99*'-'
474        for Ename in FFtable:
475            ffdata = FFtable[Ename]
476            fa = ffdata['fa']
477            fb = ffdata['fb']
478            print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
479                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
480               
481    def PrintBLtable(BLtable):
482        print '\n Neutron scattering factors:'
483        print '   Symbol   isotope       mass       b       resonant terms'
484        print 99*'-'
485        for Ename in BLtable:
486            bldata = BLtable[Ename]
487            isotope = bldata[0]
488            mass = bldata[1][0]
489            blen = bldata[1][1]
490            bres = []
491            if len(bldata[1]) > 2:
492                bres = bldata[1][2:]
493            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
494            for item in bres:
495                line += '%10.5g'%(item)
496            print line
497               
498    def PrintAtoms(General,Atoms):
499        print '\n Atoms:'
500        line = '   name    type  refine?   x         y         z    '+ \
501            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
502        if General['Type'] == 'magnetic':
503            line += '   Mx     My     Mz'
504        elif General['Type'] == 'macromolecular':
505            line = ' res no  residue  chain '+line
506        print line
507        if General['Type'] == 'nuclear':
508            print 135*'-'
509            for i,at in enumerate(Atoms):
510                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
511                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
512                if at[9] == 'I':
513                    line += '%8.4f'%(at[10])+48*' '
514                else:
515                    line += 8*' '
516                    for j in range(6):
517                        line += '%8.4f'%(at[11+j])
518                print line
519       
520    def PrintTexture(textureData):
521        topstr = '\n Spherical harmonics texture: Order:' + \
522            str(textureData['Order'])
523        if textureData['Order']:
524            print topstr+' Refine? '+str(textureData['SH Coeff'][0])
525        else:
526            print topstr
527            return
528        names = ['omega','chi','phi']
529        line = '\n'
530        for name in names:
531            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
532        print line
533        print '\n Texture coefficients:'
534        ptlbls = ' names :'
535        ptstr =  ' values:'
536        SHcoeff = textureData['SH Coeff'][1]
537        for item in SHcoeff:
538            ptlbls += '%12s'%(item)
539            ptstr += '%12.4f'%(SHcoeff[item]) 
540        print ptlbls
541        print ptstr   
542       
543    if Print: print ' Phases:'
544    phaseVary = []
545    phaseDict = {}
546    phaseConstr = {}
547    pawleyLookup = {}
548    FFtables = {}                   #scattering factors - xrays
549    BLtables = {}                   # neutrons
550    Natoms = {}
551    AtMults = {}
552    AtIA = {}
553    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
554    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
555    for name in PhaseData:
556        General = PhaseData[name]['General']
557        pId = PhaseData[name]['pId']
558        pfx = str(pId)+'::'
559        FFtable = GetFFtable(General)
560        BLtable = GetBLtable(General)
561        FFtables.update(FFtable)
562        BLtables.update(BLtable)
563        Atoms = PhaseData[name]['Atoms']
564        try:
565            PawleyRef = PhaseData[name]['Pawley ref']
566        except KeyError:
567            PawleyRef = []
568        SGData = General['SGData']
569        SGtext = G2spc.SGPrint(SGData)
570        cell = General['Cell']
571        A = G2lat.cell2A(cell[1:7])
572        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]})
573        if cell[0]:
574            phaseVary += cellVary(pfx,SGData)
575        Natoms[pfx] = 0
576        if Atoms:
577            if General['Type'] == 'nuclear':
578                Natoms[pfx] = len(Atoms)
579                for i,at in enumerate(Atoms):
580                    phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
581                        pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
582                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
583                        pfx+'AI/A:'+str(i):at[9],})
584                    if at[9] == 'I':
585                        phaseDict[pfx+'AUiso:'+str(i)] = at[10]
586                    else:
587                        phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
588                            pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
589                    if 'F' in at[2]:
590                        phaseVary.append(pfx+'Afrac:'+str(i))
591                    if 'X' in at[2]:
592                        xId,xCoef = G2spc.GetCSxinel(at[7])
593                        delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
594                        for j in range(3):
595                            if xId[j] > 0:                               
596                                phaseVary.append(delnames[j])
597                                for k in range(j):
598                                    if xId[j] == xId[k]:
599                                        G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),)) 
600                    if 'U' in at[2]:
601                        if at[9] == 'I':
602                            phaseVary.append(pfx+'AUiso:'+str(i))
603                        else:
604                            uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
605                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
606                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
607                            for j in range(6):
608                                if uId[j] > 0:                               
609                                    phaseVary.append(names[j])
610                                    for k in range(j):
611                                        if uId[j] == uId[k]:
612                                            G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
613#            elif General['Type'] == 'magnetic':
614#            elif General['Type'] == 'macromolecular':
615
616               
617            if 'SH Texture' in General:
618                textureData = General['SH Texture']
619                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
620                phaseDict[pfx+'SHorder'] = textureData['Order']
621                for name in ['omega','chi','phi']:
622                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
623                    if textureData['Sample '+name][0]:
624                        phaseVary.append(pfx+'SH '+name)
625                for name in textureData['SH Coeff'][1]:
626                    phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
627                    if textureData['SH Coeff'][0]:
628                        phaseVary.append(pfx+name)
629               
630            if Print:
631                print '\n Phase name: ',General['Name']
632                print 135*'-'
633                PrintFFtable(FFtable)
634                PrintBLtable(BLtable)
635                print ''
636                for line in SGtext: print line
637                PrintAtoms(General,Atoms)
638                print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
639                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
640                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
641                if 'SH Texture' in General:
642                    PrintTexture(textureData)
643                   
644        elif PawleyRef:
645            pawleyVary = []
646            for i,refl in enumerate(PawleyRef):
647                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
648                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
649                if refl[5]:
650                    pawleyVary.append(pfx+'PWLref:'+str(i))
651            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
652            phaseVary += pawleyVary
653               
654    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
655   
656def getVCov(varyNames,varyList,covMatrix):
657    vcov = np.zeros((len(varyNames),len(varyNames)))
658    for i1,name1 in enumerate(varyNames):
659        for i2,name2 in enumerate(varyNames):
660            try:
661                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
662            except ValueError:
663                vcov[i1][i2] = 0.0
664    return vcov
665   
666def cellFill(pfx,SGData,parmDict,sigDict): 
667    if SGData['SGLaue'] in ['-1',]:
668        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
669            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
670        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
671            sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
672    elif SGData['SGLaue'] in ['2/m',]:
673        if SGData['SGUniq'] == 'a':
674            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
675                parmDict[pfx+'A3'],0,0]
676            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
677                sigDict[pfx+'A3'],0,0]
678        elif SGData['SGUniq'] == 'b':
679            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
680                0,parmDict[pfx+'A4'],0]
681            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
682                0,sigDict[pfx+'A4'],0]
683        else:
684            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
685                0,0,parmDict[pfx+'A5']]
686            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
687                0,0,sigDict[pfx+'A5']]
688    elif SGData['SGLaue'] in ['mmm',]:
689        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
690        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
691    elif SGData['SGLaue'] in ['4/m','4/mmm']:
692        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
693        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
694    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
695        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
696            parmDict[pfx+'A0'],0,0]
697        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
698    elif SGData['SGLaue'] in ['3R', '3mR']:
699        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
700            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
701        sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
702    elif SGData['SGLaue'] in ['m3m','m3']:
703        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
704        sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
705    return A,sigA
706       
707def getCellEsd(pfx,SGData,A,covData):
708    dpr = 180./np.pi
709    rVsq = G2lat.calc_rVsq(A)
710    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
711    cell = np.array(G2lat.Gmat2cell(g))   #real cell
712    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
713    scos = cosd(cellst[3:6])
714    ssin = sind(cellst[3:6])
715    scot = scos/ssin
716    rcos = cosd(cell[3:6])
717    rsin = sind(cell[3:6])
718    rcot = rcos/rsin
719    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
720    varyList = covData['varyList']
721    covMatrix = covData['covMatrix']
722    vcov = getVCov(RMnames,varyList,covMatrix)
723    Ax = np.array(A)
724    Ax[3:] /= 2.
725    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
726        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
727    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
728    Vol = 1/np.sqrt(rVsq)
729    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
730    R123 = Ax[0]*Ax[1]*Ax[2]
731    dsasdg = np.zeros((3,6))
732    dadg = np.zeros((6,6))
733    for i0 in range(3):         #0  1   2
734        i1 = (i0+1)%3           #1  2   0
735        i2 = (i1+1)%3           #2  0   1
736        i3 = 5-i2               #3  5   4
737        i4 = 5-i1               #4  3   5
738        i5 = 5-i0               #5  4   3
739        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
740        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
741        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
742        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
743        denom = np.sqrt(denmsq)
744        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
745        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
746        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
747        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
748        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
749        dadg[i5][i5] = -Ax[i0]/denom
750    for i0 in range(3):
751        i1 = (i0+1)%3
752        i2 = (i1+1)%3
753        i3 = 5-i2
754        for ij in range(6):
755            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
756            if ij == i0:
757                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
758            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
759    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
760    var = np.diag(sigMat)
761    CS = np.where(var>0.,np.sqrt(var),0.)
762    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
763    return cellSig           
764   
765def SetPhaseData(parmDict,sigDict,Phases,covData):
766   
767    def PrintAtomsAndSig(General,Atoms,atomsSig):
768        print '\n Atoms:'
769        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
770        if General['Type'] == 'magnetic':
771            line += '   Mx     My     Mz'
772        elif General['Type'] == 'macromolecular':
773            line = ' res no  residue  chain '+line
774        print line
775        if General['Type'] == 'nuclear':
776            print 135*'-'
777            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
778                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
779            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
780            for i,at in enumerate(Atoms):
781                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
782                valstr = ' values:'
783                sigstr = ' sig   :'
784                for ind in [3,4,5,6]:
785                    sigind = str(i)+':'+str(ind)
786                    valstr += fmt[ind]%(at[ind])                   
787                    if sigind in atomsSig:
788                        sigstr += fmt[ind]%(atomsSig[sigind])
789                    else:
790                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
791                if at[9] == 'I':
792                    valstr += fmt[10]%(at[10])
793                    if str(i)+':10' in atomsSig:
794                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
795                    else:
796                        sigstr += 8*' '
797                else:
798                    valstr += 8*' '
799                    sigstr += 8*' '
800                    for ind in [11,12,13,14,15,16]:
801                        sigind = str(i)+':'+str(ind)
802                        valstr += fmt[ind]%(at[ind])
803                        if sigind in atomsSig:                       
804                            sigstr += fmt[ind]%(atomsSig[sigind])
805                        else:
806                            sigstr += 8*' '
807                print name
808                print valstr
809                print sigstr
810               
811    def PrintSHtextureAndSig(textureData,SHtextureSig):
812        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
813        names = ['omega','chi','phi']
814        namstr = '  names :'
815        ptstr =  '  values:'
816        sigstr = '  esds  :'
817        for name in names:
818            namstr += '%12s'%(name)
819            ptstr += '%12.3f'%(textureData['Sample '+name][1])
820            if 'Sample '+name in SHtextureSig:
821                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
822            else:
823                sigstr += 12*' '
824        print namstr
825        print ptstr
826        print sigstr
827        print '\n Texture coefficients:'
828        namstr = '  names :'
829        ptstr =  '  values:'
830        sigstr = '  esds  :'
831        SHcoeff = textureData['SH Coeff'][1]
832        for name in SHcoeff:
833            namstr += '%12s'%(name)
834            ptstr += '%12.3f'%(SHcoeff[name])
835            if name in SHtextureSig:
836                sigstr += '%12.3f'%(SHtextureSig[name])
837            else:
838                sigstr += 12*' '
839        print namstr
840        print ptstr
841        print sigstr
842       
843           
844    print '\n Phases:'
845    for phase in Phases:
846        print ' Result for phase: ',phase
847        Phase = Phases[phase]
848        General = Phase['General']
849        SGData = General['SGData']
850        Atoms = Phase['Atoms']
851        cell = General['Cell']
852        pId = Phase['pId']
853        pfx = str(pId)+'::'
854        if cell[0]:
855            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
856            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
857            print ' Reciprocal metric tensor: '
858            ptfmt = "%15.9f"
859            names = ['A11','A22','A33','A12','A13','A23']
860            namstr = '  names :'
861            ptstr =  '  values:'
862            sigstr = '  esds  :'
863            for name,a,siga in zip(names,A,sigA):
864                namstr += '%15s'%(name)
865                ptstr += ptfmt%(a)
866                if siga:
867                    sigstr += ptfmt%(siga)
868                else:
869                    sigstr += 15*' '
870            print namstr
871            print ptstr
872            print sigstr
873            cell[1:7] = G2lat.A2cell(A)
874            cell[7] = G2lat.calc_V(A)
875            print ' New unit cell:'
876            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
877            names = ['a','b','c','alpha','beta','gamma','Volume']
878            namstr = '  names :'
879            ptstr =  '  values:'
880            sigstr = '  esds  :'
881            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
882                namstr += '%12s'%(name)
883                ptstr += fmt%(a)
884                if siga:
885                    sigstr += fmt%(siga)
886                else:
887                    sigstr += 12*' '
888            print namstr
889            print ptstr
890            print sigstr
891           
892        if 'Pawley' in Phase['General']['Type']:
893            pawleyRef = Phase['Pawley ref']
894            for i,refl in enumerate(pawleyRef):
895                key = pfx+'PWLref:'+str(i)
896                refl[6] = abs(parmDict[key])        #suppress negative Fsq
897                if key in sigDict:
898                    refl[7] = sigDict[key]
899                else:
900                    refl[7] = 0
901        else:
902            atomsSig = {}
903            if General['Type'] == 'nuclear':
904                for i,at in enumerate(Atoms):
905                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
906                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
907                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
908                    for ind in [3,4,5,6]:
909                        at[ind] = parmDict[names[ind]]
910                        if ind in [3,4,5]:
911                            name = names[ind].replace('A','dA')
912                        else:
913                            name = names[ind]
914                        if name in sigDict:
915                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
916                    if at[9] == 'I':
917                        at[10] = parmDict[names[10]]
918                        if names[10] in sigDict:
919                            atomsSig[str(i)+':10'] = sigDict[names[10]]
920                    else:
921                        for ind in [11,12,13,14,15,16]:
922                            at[ind] = parmDict[names[ind]]
923                            if names[ind] in sigDict:
924                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
925            PrintAtomsAndSig(General,Atoms,atomsSig)
926       
927        if 'SH Texture' in General:
928            textureData = General['SH Texture']   
929            if textureData['Order']:
930                SHtextureSig = {}
931                for name in ['omega','chi','phi']:
932                    aname = pfx+'SH '+name
933                    textureData['Sample '+name][1] = parmDict[aname]
934                    if aname in sigDict:
935                        SHtextureSig['Sample '+name] = sigDict[aname]
936                for name in textureData['SH Coeff'][1]:
937                    aname = pfx+name
938                    textureData['SH Coeff'][1][name] = parmDict[aname]
939                    if aname in sigDict:
940                        SHtextureSig[name] = sigDict[aname]
941                PrintSHtextureAndSig(textureData,SHtextureSig)
942
943def GetHistogramPhaseData(Phases,Histograms,Print=True):
944   
945    def PrintSize(hapData):
946        if hapData[0] in ['isotropic','uniaxial']:
947            line = '\n Size model    : %9s'%(hapData[0])
948            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
949            if hapData[0] == 'uniaxial':
950                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
951            print line
952        else:
953            print '\n Size model    : %s'%(hapData[0])
954            Snames = ['S11','S22','S33','S12','S13','S23']
955            ptlbls = ' names :'
956            ptstr =  ' values:'
957            varstr = ' refine:'
958            for i,name in enumerate(Snames):
959                ptlbls += '%12s' % (name)
960                ptstr += '%12.6f' % (hapData[4][i])
961                varstr += '%12s' % (str(hapData[5][i]))
962            print ptlbls
963            print ptstr
964            print varstr
965       
966    def PrintMuStrain(hapData,SGData):
967        if hapData[0] in ['isotropic','uniaxial']:
968            line = '\n Mustrain model: %9s'%(hapData[0])
969            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
970            if hapData[0] == 'uniaxial':
971                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
972            print line
973        else:
974            print '\n Mustrain model: %s'%(hapData[0])
975            Snames = G2spc.MustrainNames(SGData)
976            ptlbls = ' names :'
977            ptstr =  ' values:'
978            varstr = ' refine:'
979            for i,name in enumerate(Snames):
980                ptlbls += '%12s' % (name)
981                ptstr += '%12.6f' % (hapData[4][i])
982                varstr += '%12s' % (str(hapData[5][i]))
983            print ptlbls
984            print ptstr
985            print varstr
986
987    def PrintHStrain(hapData,SGData):
988        print '\n Hydrostatic/elastic strain: '
989        Hsnames = G2spc.HStrainNames(SGData)
990        ptlbls = ' names :'
991        ptstr =  ' values:'
992        varstr = ' refine:'
993        for i,name in enumerate(Hsnames):
994            ptlbls += '%12s' % (name)
995            ptstr += '%12.6f' % (hapData[0][i])
996            varstr += '%12s' % (str(hapData[1][i]))
997        print ptlbls
998        print ptstr
999        print varstr
1000
1001    def PrintSHPO(hapData):
1002        print '\n Spherical harmonics preferred orientation: Order:' + \
1003            str(hapData[4])+' Refine? '+str(hapData[2])
1004        ptlbls = ' names :'
1005        ptstr =  ' values:'
1006        for item in hapData[5]:
1007            ptlbls += '%12s'%(item)
1008            ptstr += '%12.3f'%(hapData[5][item]) 
1009        print ptlbls
1010        print ptstr
1011   
1012    hapDict = {}
1013    hapVary = []
1014    controlDict = {}
1015    poType = {}
1016    poAxes = {}
1017    spAxes = {}
1018    spType = {}
1019   
1020    for phase in Phases:
1021        HistoPhase = Phases[phase]['Histograms']
1022        SGData = Phases[phase]['General']['SGData']
1023        cell = Phases[phase]['General']['Cell'][1:7]
1024        A = G2lat.cell2A(cell)
1025        pId = Phases[phase]['pId']
1026        for histogram in HistoPhase:
1027            try:
1028                Histogram = Histograms[histogram]
1029            except KeyError:                       
1030                #skip if histogram not included e.g. in a sequential refinement
1031                continue
1032            hapData = HistoPhase[histogram]
1033            hId = Histogram['hId']
1034            limits = Histogram['Limits'][1]
1035            inst = Histogram['Instrument Parameters']
1036            inst = dict(zip(inst[3],inst[1]))
1037            Zero = inst['Zero']
1038            if 'C' in inst['Type']:
1039                try:
1040                    wave = inst['Lam']
1041                except KeyError:
1042                    wave = inst['Lam1']
1043                dmin = wave/(2.0*sind(limits[1]/2.0))
1044            pfx = str(pId)+':'+str(hId)+':'
1045            for item in ['Scale','Extinction']:
1046                hapDict[pfx+item] = hapData[item][0]
1047                if hapData[item][1]:
1048                    hapVary.append(pfx+item)
1049            names = G2spc.HStrainNames(SGData)
1050            for i,name in enumerate(names):
1051                hapDict[pfx+name] = hapData['HStrain'][0][i]
1052                if hapData['HStrain'][1][i]:
1053                    hapVary.append(pfx+name)
1054            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1055            if hapData['Pref.Ori.'][0] == 'MD':
1056                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1057                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1058                if hapData['Pref.Ori.'][2]:
1059                    hapVary.append(pfx+'MD')
1060            else:                           #'SH' spherical harmonics
1061                controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1062                controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1063                for item in hapData['Pref.Ori.'][5]:
1064                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1065                    if hapData['Pref.Ori.'][2]:
1066                        hapVary.append(pfx+item)
1067            for item in ['Mustrain','Size']:
1068                controlDict[pfx+item+'Type'] = hapData[item][0]
1069                if hapData[item][0] in ['isotropic','uniaxial']:
1070                    hapDict[pfx+item+':0'] = hapData[item][1][0]
1071                    if hapData[item][2][0]:
1072                        hapVary.append(pfx+item+':0')
1073                    if hapData[item][0] == 'uniaxial':
1074                        controlDict[pfx+item+'Axis'] = hapData[item][3]
1075                        hapDict[pfx+item+':1'] = hapData[item][1][1]
1076                        if hapData[item][2][1]:
1077                            hapVary.append(pfx+item+':1')
1078                else:       #generalized for mustrain or ellipsoidal for size
1079                    if item == 'Mustrain':
1080                        names = G2spc.MustrainNames(SGData)
1081                        pwrs = []
1082                        for name in names:
1083                            h,k,l = name[1:]
1084                            pwrs.append([int(h),int(k),int(l)])
1085                        controlDict[pfx+'MuPwrs'] = pwrs
1086                    for i in range(len(hapData[item][4])):
1087                        sfx = ':'+str(i)
1088                        hapDict[pfx+item+sfx] = hapData[item][4][i]
1089                        if hapData[item][5][i]:
1090                            hapVary.append(pfx+item+sfx)
1091                           
1092            if Print: 
1093                print '\n Phase: ',phase,' in histogram: ',histogram
1094                print 135*'-'
1095                print ' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1096                print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1097                if hapData['Pref.Ori.'][0] == 'MD':
1098                    Ax = hapData['Pref.Ori.'][3]
1099                    print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1100                        ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1101                else: #'SH' for spherical harmonics
1102                    PrintSHPO(hapData['Pref.Ori.'])
1103                PrintSize(hapData['Size'])
1104                PrintMuStrain(hapData['Mustrain'],SGData)
1105                PrintHStrain(hapData['HStrain'],SGData)
1106            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1107            refList = []
1108            for h,k,l,d in HKLd:
1109                ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1110                if ext:
1111                    continue
1112                if 'C' in inst['Type']:
1113                    pos = 2.0*asind(wave/(2.0*d))
1114                    if limits[0] < pos < limits[1]:
1115                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0])
1116                else:
1117                    raise ValueError 
1118            Histogram['Reflection Lists'][phase] = refList
1119    return hapVary,hapDict,controlDict
1120   
1121def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
1122   
1123    def PrintSizeAndSig(hapData,sizeSig):
1124        line = '\n Size model:     %9s'%(hapData[0])
1125        if hapData[0] in ['isotropic','uniaxial']:
1126            line += ' equatorial:%12.3f'%(hapData[1][0])
1127            if sizeSig[0][0]:
1128                line += ', sig: %8.3f'%(sizeSig[0][0])
1129            if hapData[0] == 'uniaxial':
1130                line += ' axial:%12.3f'%(hapData[1][1])
1131                if sizeSig[0][1]:
1132                    line += ', sig: %8.3f'%(sizeSig[0][1])
1133            print line
1134        else:
1135            print line
1136            Snames = ['S11','S22','S33','S12','S13','S23']
1137            ptlbls = ' name  :'
1138            ptstr =  ' value :'
1139            sigstr = ' sig   :'
1140            for i,name in enumerate(Snames):
1141                ptlbls += '%12s' % (name)
1142                ptstr += '%12.6f' % (hapData[4][i])
1143                if sizeSig[1][i]:
1144                    sigstr += '%12.6f' % (sizeSig[1][i])
1145                else:
1146                    sigstr += 12*' '
1147            print ptlbls
1148            print ptstr
1149            print sigstr
1150       
1151    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1152        line = '\n Mustrain model: %9s'%(hapData[0])
1153        if hapData[0] in ['isotropic','uniaxial']:
1154            line += ' equatorial:%12.1f'%(hapData[1][0])
1155            if mustrainSig[0][0]:
1156                line += ', sig: %8.1f'%(mustrainSig[0][0])
1157            if hapData[0] == 'uniaxial':
1158                line += ' axial:%12.1f'%(hapData[1][1])
1159                if mustrainSig[0][1]:
1160                     line += ', sig: %8.1f'%(mustrainSig[0][1])
1161            print line
1162        else:
1163            print line
1164            Snames = G2spc.MustrainNames(SGData)
1165            ptlbls = ' name  :'
1166            ptstr =  ' value :'
1167            sigstr = ' sig   :'
1168            for i,name in enumerate(Snames):
1169                ptlbls += '%12s' % (name)
1170                ptstr += '%12.6f' % (hapData[4][i])
1171                if mustrainSig[1][i]:
1172                    sigstr += '%12.6f' % (mustrainSig[1][i])
1173                else:
1174                    sigstr += 12*' '
1175            print ptlbls
1176            print ptstr
1177            print sigstr
1178           
1179    def PrintHStrainAndSig(hapData,strainSig,SGData):
1180        print '\n Hydrostatic/elastic strain: '
1181        Hsnames = G2spc.HStrainNames(SGData)
1182        ptlbls = ' name  :'
1183        ptstr =  ' value :'
1184        sigstr = ' sig   :'
1185        for i,name in enumerate(Hsnames):
1186            ptlbls += '%12s' % (name)
1187            ptstr += '%12.6g' % (hapData[0][i])
1188            if name in strainSig:
1189                sigstr += '%12.6g' % (strainSig[name])
1190            else:
1191                sigstr += 12*' '
1192        print ptlbls
1193        print ptstr
1194        print sigstr
1195       
1196    def PrintSHPOAndSig(hapData,POsig):
1197        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1198        ptlbls = ' names :'
1199        ptstr =  ' values:'
1200        sigstr = ' sig   :'
1201        for item in hapData[5]:
1202            ptlbls += '%12s'%(item)
1203            ptstr += '%12.3f'%(hapData[5][item])
1204            if item in POsig:
1205                sigstr += '%12.3f'%(POsig[item])
1206            else:
1207                sigstr += 12*' ' 
1208        print ptlbls
1209        print ptstr
1210        print sigstr
1211   
1212    for phase in Phases:
1213        HistoPhase = Phases[phase]['Histograms']
1214        SGData = Phases[phase]['General']['SGData']
1215        pId = Phases[phase]['pId']
1216        for histogram in HistoPhase:
1217            try:
1218                Histogram = Histograms[histogram]
1219            except KeyError:                       
1220                #skip if histogram not included e.g. in a sequential refinement
1221                continue
1222            print '\n Phase: ',phase,' in histogram: ',histogram
1223            print 130*'-'
1224            hapData = HistoPhase[histogram]
1225            hId = Histogram['hId']
1226            pfx = str(pId)+':'+str(hId)+':'
1227            print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
1228                %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1229           
1230            PhFrExtPOSig = {}
1231            for item in ['Scale','Extinction']:
1232                hapData[item][0] = parmDict[pfx+item]
1233                if pfx+item in sigDict:
1234                    PhFrExtPOSig[item] = sigDict[pfx+item]
1235            if hapData['Pref.Ori.'][0] == 'MD':
1236                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1237                if pfx+'MD' in sigDict:
1238                    PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
1239            else:                           #'SH' spherical harmonics
1240                for item in hapData['Pref.Ori.'][5]:
1241                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1242                    if pfx+item in sigDict:
1243                        PhFrExtPOSig[item] = sigDict[pfx+item]
1244            if Print:
1245                if 'Scale' in PhFrExtPOSig:
1246                    print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
1247                if 'Extinction' in PhFrExtPOSig:
1248                    print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
1249                if hapData['Pref.Ori.'][0] == 'MD':
1250                    if 'MD' in PhFrExtPOSig:
1251                        print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
1252                else:
1253                    PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
1254            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1255                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
1256                'HStrain':{}}                 
1257            for item in ['Mustrain','Size']:
1258                if hapData[item][0] in ['isotropic','uniaxial']:                   
1259                    hapData[item][1][0] = parmDict[pfx+item+':0']
1260                    if item == 'Size':
1261                        hapData[item][1][0] = min(10.,max(0.01,hapData[item][1][0]))
1262                    if pfx+item+':0' in sigDict: 
1263                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':0']
1264                    if hapData[item][0] == 'uniaxial':
1265                        hapData[item][1][1] = parmDict[pfx+item+':1']
1266                        if item == 'Size':
1267                            hapData[item][1][1] = min(10.,max(0.01,hapData[item][1][1]))                       
1268                        if pfx+item+':1' in sigDict:
1269                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':1']
1270                else:       #generalized for mustrain or ellipsoidal for size
1271                    for i in range(len(hapData[item][4])):
1272                        sfx = ':'+str(i)
1273                        hapData[item][4][i] = parmDict[pfx+item+sfx]
1274                        if pfx+item+sfx in sigDict:
1275                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
1276            names = G2spc.HStrainNames(SGData)
1277            for i,name in enumerate(names):
1278                hapData['HStrain'][0][i] = parmDict[pfx+name]
1279                if pfx+name in sigDict:
1280                    SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
1281            if Print:
1282                PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
1283                PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
1284                PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
1285   
1286def GetHistogramData(Histograms,Print=True):
1287   
1288    def GetBackgroundParms(hId,Background):
1289        bakType,bakFlag = Background[:2]
1290        backVals = Background[3:]
1291        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1292        if bakFlag:                                 #returns backNames as varyList = backNames
1293            return bakType,dict(zip(backNames,backVals)),backNames
1294        else:                                       #no background varied; varyList = []
1295            return bakType,dict(zip(backNames,backVals)),[]
1296       
1297    def GetInstParms(hId,Inst):
1298        insVals,insFlags,insNames = Inst[1:4]
1299        dataType = insVals[0]
1300        instDict = {}
1301        insVary = []
1302        pfx = ':'+str(hId)+':'
1303        for i,flag in enumerate(insFlags):
1304            insName = pfx+insNames[i]
1305            instDict[insName] = insVals[i]
1306            if flag:
1307                insVary.append(insName)
1308        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.01)
1309        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.01)
1310        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1311        return dataType,instDict,insVary
1312       
1313    def GetSampleParms(hId,Sample):
1314        sampVary = []
1315        hfx = ':'+str(hId)+':'       
1316        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1317            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1318        Type = Sample['Type']
1319        if 'Bragg' in Type:             #Bragg-Brentano
1320            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
1321                sampDict[hfx+item] = Sample[item][0]
1322                if Sample[item][1]:
1323                    sampVary.append(hfx+item)
1324        elif 'Debye' in Type:        #Debye-Scherrer
1325            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1326                sampDict[hfx+item] = Sample[item][0]
1327                if Sample[item][1]:
1328                    sampVary.append(hfx+item)
1329        return Type,sampDict,sampVary
1330       
1331    def PrintBackground(Background):
1332        print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
1333        line = ' Coefficients: '
1334        for i,back in enumerate(Background[3:]):
1335            line += '%10.3f'%(back)
1336            if i and not i%10:
1337                line += '\n'+15*' '
1338        print line
1339       
1340    def PrintInstParms(Inst):
1341        print '\n Instrument Parameters:'
1342        ptlbls = ' name  :'
1343        ptstr =  ' value :'
1344        varstr = ' refine:'
1345        instNames = Inst[3][1:]
1346        for i,name in enumerate(instNames):
1347            ptlbls += '%12s' % (name)
1348            ptstr += '%12.6f' % (Inst[1][i+1])
1349            if name in ['Lam1','Lam2','Azimuth']:
1350                varstr += 12*' '
1351            else:
1352                varstr += '%12s' % (str(bool(Inst[2][i+1])))
1353        print ptlbls
1354        print ptstr
1355        print varstr
1356       
1357    def PrintSampleParms(Sample):
1358        print '\n Sample Parameters:'
1359        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1360            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
1361        ptlbls = ' name  :'
1362        ptstr =  ' value :'
1363        varstr = ' refine:'
1364        if 'Bragg' in Sample['Type']:
1365            for item in ['Scale','Shift','Transparency']:
1366                ptlbls += '%14s'%(item)
1367                ptstr += '%14.4f'%(Sample[item][0])
1368                varstr += '%14s'%(str(bool(Sample[item][1])))
1369           
1370        elif 'Debye' in Type:        #Debye-Scherrer
1371            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1372                ptlbls += '%14s'%(item)
1373                ptstr += '%14.4f'%(Sample[item][0])
1374                varstr += '%14s'%(str(bool(Sample[item][1])))
1375
1376        print ptlbls
1377        print ptstr
1378        print varstr
1379       
1380
1381    histDict = {}
1382    histVary = []
1383    controlDict = {}
1384    for histogram in Histograms:
1385        Histogram = Histograms[histogram]
1386        hId = Histogram['hId']
1387        pfx = ':'+str(hId)+':'
1388        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1389       
1390        Background = Histogram['Background'][0]
1391        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1392        controlDict[pfx+'bakType'] = Type
1393        histDict.update(bakDict)
1394        histVary += bakVary
1395       
1396        Inst = Histogram['Instrument Parameters']
1397        Type,instDict,insVary = GetInstParms(hId,Inst)
1398        controlDict[pfx+'histType'] = Type
1399        if pfx+'Lam1' in instDict:
1400            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1401        else:
1402            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
1403        histDict.update(instDict)
1404        histVary += insVary
1405       
1406        Sample = Histogram['Sample Parameters']
1407        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1408        controlDict[pfx+'instType'] = Type
1409        histDict.update(sampDict)
1410        histVary += sampVary
1411
1412        if Print: 
1413            print '\n Histogram: ',histogram,' histogram Id: ',hId
1414            print 135*'-'
1415            Units = {'C':' deg','T':' msec'}
1416            units = Units[controlDict[pfx+'histType'][2]]
1417            Limits = controlDict[pfx+'Limits']
1418            print ' Instrument type: ',Sample['Type']
1419            print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
1420            PrintSampleParms(Sample)
1421            PrintInstParms(Inst)
1422            PrintBackground(Background)
1423       
1424    return histVary,histDict,controlDict
1425   
1426def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
1427   
1428    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1429        lenBack = len(Background[3:])
1430        backSig = [0 for i in range(lenBack)]
1431        for i in range(lenBack):
1432            Background[3+i] = parmDict[pfx+'Back:'+str(i)]
1433            if pfx+'Back:'+str(i) in sigDict:
1434                backSig[i] = sigDict[pfx+'Back:'+str(i)]
1435        return backSig
1436       
1437    def SetInstParms(pfx,Inst,parmDict,sigDict):
1438        insVals,insFlags,insNames = Inst[1:4]
1439        instSig = [0 for i in range(len(insVals))]
1440        for i,flag in enumerate(insFlags):
1441            insName = pfx+insNames[i]
1442            insVals[i] = parmDict[insName]
1443            if insName in sigDict:
1444                instSig[i] = sigDict[insName]
1445        return instSig
1446       
1447    def SetSampleParms(pfx,Sample,parmDict,sigDict):
1448        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
1449            sampSig = [0 for i in range(3)]
1450            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
1451                Sample[item][0] = parmDict[pfx+item]
1452                if pfx+item in sigDict:
1453                    sampSig[i] = sigDict[pfx+item]
1454        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1455            sampSig = [0 for i in range(4)]
1456            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1457                Sample[item][0] = parmDict[pfx+item]
1458                if pfx+item in sigDict:
1459                    sampSig[i] = sigDict[pfx+item]
1460        return sampSig
1461       
1462    def PrintBackgroundSig(Background,backSig):
1463        print '\n Background function: ',Background[0]
1464        valstr = ' value : '
1465        sigstr = ' sig   : '
1466        for i,back in enumerate(Background[3:]):
1467            valstr += '%10.4f'%(back)
1468            if Background[1]:
1469                sigstr += '%10.4f'%(backSig[i])
1470            else:
1471                sigstr += 10*' '
1472        print valstr
1473        print sigstr
1474       
1475    def PrintInstParmsSig(Inst,instSig):
1476        print '\n Instrument Parameters:'
1477        ptlbls = ' names :'
1478        ptstr =  ' value :'
1479        sigstr = ' sig   :'
1480        instNames = Inst[3][1:]
1481        for i,name in enumerate(instNames):
1482            ptlbls += '%12s' % (name)
1483            ptstr += '%12.6f' % (Inst[1][i+1])
1484            if instSig[i+1]:
1485                sigstr += '%12.6f' % (instSig[i+1])
1486            else:
1487                sigstr += 12*' '
1488        print ptlbls
1489        print ptstr
1490        print sigstr
1491       
1492    def PrintSampleParmsSig(Sample,sampleSig):
1493        print '\n Sample Parameters:'
1494        ptlbls = ' names :'
1495        ptstr =  ' values:'
1496        sigstr = ' sig   :'
1497        if 'Bragg' in Sample['Type']:
1498            for i,item in enumerate(['Scale','Shift','Transparency']):
1499                ptlbls += '%14s'%(item)
1500                ptstr += '%14.4f'%(Sample[item][0])
1501                if sampleSig[i]:
1502                    sigstr += '%14.4f'%(sampleSig[i])
1503                else:
1504                    sigstr += 14*' '
1505           
1506        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1507            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1508                ptlbls += '%14s'%(item)
1509                ptstr += '%14.4f'%(Sample[item][0])
1510                if sampleSig[i]:
1511                    sigstr += '%14.4f'%(sampleSig[i])
1512                else:
1513                    sigstr += 14*' '
1514
1515        print ptlbls
1516        print ptstr
1517        print sigstr
1518       
1519    for histogram in Histograms:
1520        if 'PWDR' in histogram:
1521            Histogram = Histograms[histogram]
1522            hId = Histogram['hId']
1523            pfx = ':'+str(hId)+':'
1524            Background = Histogram['Background'][0]
1525            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1526           
1527            Inst = Histogram['Instrument Parameters']
1528            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1529       
1530            Sample = Histogram['Sample Parameters']
1531            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1532
1533            print '\n Histogram: ',histogram,' histogram Id: ',hId
1534            print 135*'-'
1535            print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
1536            if Print:
1537                print ' Instrument type: ',Sample['Type']
1538                PrintSampleParmsSig(Sample,sampSig)
1539                PrintInstParmsSig(Inst,instSig)
1540                PrintBackgroundSig(Background,backSig)
1541
1542def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict):
1543    Natoms = calcControls['Natoms'][pfx]
1544    Tdata = Natoms*[' ',]
1545    Mdata = np.zeros(Natoms)
1546    IAdata = Natoms*[' ',]
1547    Fdata = np.zeros(Natoms)
1548    FFdata = []
1549    BLdata = []
1550    Xdata = np.zeros((3,Natoms))
1551    dXdata = np.zeros((3,Natoms))
1552    Uisodata = np.zeros(Natoms)
1553    Uijdata = np.zeros((6,Natoms))
1554    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1555        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1556        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1557        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1558        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1559    for iatm in range(Natoms):
1560        for key in keys:
1561            parm = pfx+key+str(iatm)
1562            if parm in parmDict:
1563                keys[key][iatm] = parmDict[parm]
1564        FFdata.append(FFtables[Tdata[iatm]])
1565        BLdata.append(BLtables[Tdata[iatm]][1])
1566    return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1567   
1568def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1569    ''' Compute structure factors for all h,k,l for phase
1570    input:
1571        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1572        G:      reciprocal metric tensor
1573        pfx:    phase id string
1574        SGData: space group info. dictionary output from SpcGroup
1575        calcControls:
1576        ParmDict:
1577    puts result F^2 in each ref[8] in refList
1578    '''       
1579    twopi = 2.0*np.pi
1580    twopisq = 2.0*np.pi**2
1581    ast = np.sqrt(np.diag(G))
1582    Mast = twopisq*np.multiply.outer(ast,ast)
1583    FFtables = calcControls['FFtables']
1584    BLtables = calcControls['BLtables']
1585    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
1586    if 'N' in parmDict[hfx+'Type']:
1587        FP,FPP = G2el.BlenRes(BLdata,parmDict[hfx+'Lam'])
1588    else:
1589        FP = np.array([El[hfx+'FP'] for El in FFdata])
1590        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1591    maxPos = len(SGData['SGOps'])
1592    Uij = np.array(G2lat.U6toUij(Uijdata))
1593    bij = Mast*Uij.T
1594    for refl in refList:
1595        fbs = np.array([0,0])
1596        H = refl[:3]
1597        SQ = 1./(2.*refl[4])**2
1598        if 'N' in parmDict[hfx+'Type']:
1599            FF = np.array([El[1] for El in BLdata])
1600        else:       #'X'
1601            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1602        SQfactor = 4.0*SQ*twopisq
1603        Uniq = refl[11]
1604        phi = refl[12]
1605        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1606        sinp = np.sin(phase)
1607        cosp = np.cos(phase)
1608        occ = Mdata*Fdata/len(Uniq)
1609        biso = -SQfactor*Uisodata
1610        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1611        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1612        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1613        Tcorr = Tiso*Tuij
1614        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
1615        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
1616        if not SGData['SGInv']:
1617            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
1618            fbs = np.sum(np.sum(fb,axis=1),axis=1)
1619        fasq = fas**2
1620        fbsq = fbs**2        #imaginary
1621        refl[9] = np.sum(fasq)+np.sum(fbsq)
1622        refl[10] = atan2d(fbs[0],fas[0])
1623    return refList
1624   
1625def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1626    twopi = 2.0*np.pi
1627    twopisq = 2.0*np.pi**2
1628    ast = np.sqrt(np.diag(G))
1629    Mast = twopisq*np.multiply.outer(ast,ast)
1630    FFtables = calcControls['FFtables']
1631    BLtables = calcControls['BLtables']
1632    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
1633    if 'N' in parmDict[hfx+'Type']:
1634        FP = 0.
1635        FPP = 0.
1636    else:
1637        FP = np.array([El[hfx+'FP'] for El in FFdata])
1638        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1639    maxPos = len(SGData['SGOps'])       
1640    Uij = np.array(G2lat.U6toUij(Uijdata))
1641    bij = Mast*Uij.T
1642    dFdvDict = {}
1643    dFdfr = np.zeros((len(refList),len(Mdata)))
1644    dFdx = np.zeros((len(refList),len(Mdata),3))
1645    dFdui = np.zeros((len(refList),len(Mdata)))
1646    dFdua = np.zeros((len(refList),len(Mdata),6))
1647    for iref,refl in enumerate(refList):
1648        H = np.array(refl[:3])
1649        SQ = 1./(2.*refl[4])**2          # or (sin(theta)/lambda)**2
1650        if 'N' in parmDict[hfx+'Type']:
1651            FF = np.array([El[1] for El in BLdata])
1652        else:       #'X'
1653            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1654        SQfactor = 8.0*SQ*np.pi**2
1655        Uniq = refl[11]
1656        phi = refl[12]
1657        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
1658        sinp = np.sin(phase)
1659        cosp = np.cos(phase)
1660        occ = Mdata*Fdata/len(Uniq)
1661        biso = -SQfactor*Uisodata
1662        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1663#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1664        HbH = -np.inner(H,np.inner(bij,H))
1665        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1666        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1667        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1668        Tcorr = Tiso*Tuij
1669        fot = (FF+FP)*occ*Tcorr
1670        fotp = FPP*occ*Tcorr
1671        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1672        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1673       
1674        fas = np.sum(np.sum(fa,axis=1),axis=1)
1675        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1676        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1677        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1678        #sum below is over Uniq
1679        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
1680        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1681        dfadui = np.sum(-SQfactor*fa,axis=2)
1682        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1683        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
1684        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1685        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1686        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1687        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1688        if not SGData['SGInv']:
1689            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
1690            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
1691            dfbdui = np.sum(-SQfactor*fb,axis=2)
1692            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1693            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1694            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1695            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1696            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1697        #loop over atoms - each dict entry is list of derivatives for all the reflections
1698    for i in range(len(Mdata)):     
1699        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1700        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1701        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1702        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1703        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1704        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1705        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1706        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1707        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1708        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1709        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1710    return dFdvDict
1711       
1712def Dict2Values(parmdict, varylist):
1713    '''Use before call to leastsq to setup list of values for the parameters
1714    in parmdict, as selected by key in varylist'''
1715    return [parmdict[key] for key in varylist] 
1716   
1717def Values2Dict(parmdict, varylist, values):
1718    ''' Use after call to leastsq to update the parameter dictionary with
1719    values corresponding to keys in varylist'''
1720    parmdict.update(zip(varylist,values))
1721   
1722def GetNewCellParms(parmDict,varyList):
1723    newCellDict = {}
1724    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
1725    Adict = dict(zip(['A'+str(i) for i in range(6)],['D11','D22','D33','D12','D13','D23']))
1726    for item in varyList:
1727        keys = item.split(':')
1728        if keys[2] in Ddict:
1729            key = keys[0]+':'+keys[1]+':'+Ddict[keys[2]]
1730            newCellDict[key] = parmDict[key]+parmDict[item]
1731        elif keys[2] in Adict:
1732            newCellDict[key] = parmDict[item]
1733           
1734   
1735   
1736    return newCellDict
1737   
1738def ApplyXYZshifts(parmDict,varyList):
1739    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1740        input:
1741            parmDict - parameter dictionary
1742            varyList - list of variables
1743        returns:
1744            newAtomDict - dictitemionary of new atomic coordinate names & values;
1745                key is parameter shift name
1746    '''
1747    newAtomDict = {}
1748    for item in parmDict:
1749        if 'dA' in item:
1750            parm = ''.join(item.split('d'))
1751            parmDict[parm] += parmDict[item]
1752            newAtomDict[item] = [parm,parmDict[parm]]
1753    return newAtomDict
1754   
1755def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1756    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1757    odfCor = 1.0
1758    H = refl[:3]
1759    cell = G2lat.Gmat2cell(g)
1760    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1761    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1762    phi,beta = G2lat.CrsAng(H,cell,SGData)
1763    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1764    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1765    for item in SHnames:
1766        L,M,N = eval(item.strip('C'))
1767        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1768        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1769        Lnorm = G2lat.Lnorm(L)
1770        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1771    return odfCor
1772   
1773def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1774    FORPI = 12.5663706143592
1775    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1776    odfCor = 1.0
1777    dFdODF = {}
1778    dFdSA = [0,0,0]
1779    H = refl[:3]
1780    cell = G2lat.Gmat2cell(g)
1781    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1782    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1783    phi,beta = G2lat.CrsAng(H,cell,SGData)
1784    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
1785    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1786    for item in SHnames:
1787        L,M,N = eval(item.strip('C'))
1788        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1789        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1790        Lnorm = G2lat.Lnorm(L)
1791        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1792        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1793        for i in range(3):
1794            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1795    return odfCor,dFdODF,dFdSA
1796   
1797def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1798    odfCor = 1.0
1799    H = refl[:3]
1800    cell = G2lat.Gmat2cell(g)
1801    Sangl = [0.,0.,0.]
1802    if 'Bragg' in calcControls[hfx+'instType']:
1803        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1804        IFCoup = True
1805    else:
1806        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1807        IFCoup = False
1808    phi,beta = G2lat.CrsAng(H,cell,SGData)
1809    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1810    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1811    for item in SHnames:
1812        L,N = eval(item.strip('C'))
1813        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1814        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1815    return odfCor
1816   
1817def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1818    FORPI = 12.5663706143592
1819    odfCor = 1.0
1820    dFdODF = {}
1821    H = refl[:3]
1822    cell = G2lat.Gmat2cell(g)
1823    Sangl = [0.,0.,0.]
1824    if 'Bragg' in calcControls[hfx+'instType']:
1825        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1826        IFCoup = True
1827    else:
1828        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1829        IFCoup = False
1830    phi,beta = G2lat.CrsAng(H,cell,SGData)
1831    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1832    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1833    for item in SHnames:
1834        L,N = eval(item.strip('C'))
1835        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1836        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1837        dFdODF[phfx+item] = Kcsl*Lnorm
1838    return odfCor,dFdODF
1839   
1840def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1841    if calcControls[phfx+'poType'] == 'MD':
1842        MD = parmDict[phfx+'MD']
1843        MDAxis = calcControls[phfx+'MDAxis']
1844        sumMD = 0
1845        for H in refl[11]:           
1846            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1847            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1848            sumMD += A**3
1849        POcorr = sumMD/len(refl[11])
1850    else:   #spherical harmonics
1851        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1852    return POcorr
1853   
1854def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1855    POderv = {}
1856    if calcControls[phfx+'poType'] == 'MD':
1857        MD = parmDict[phfx+'MD']
1858        MDAxis = calcControls[phfx+'MDAxis']
1859        sumMD = 0
1860        sumdMD = 0
1861        for H in refl[11]:           
1862            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1863            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1864            sumMD += A**3
1865            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1866        POcorr = sumMD/len(refl[11])
1867        POderv[phfx+'MD'] = sumdMD/len(refl[11])
1868    else:   #spherical harmonics
1869        POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1870    return POcorr,POderv
1871   
1872def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1873    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1874    if 'X' in parmDict[hfx+'Type']:
1875        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1876    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1877    if pfx+'SHorder' in parmDict:
1878        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1879    refl[13] = Icorr       
1880   
1881def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1882    dIdsh = 1./parmDict[hfx+'Scale']
1883    dIdsp = 1./parmDict[phfx+'Scale']
1884    if 'X' in parmDict[hfx+'Type']:
1885        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1886        dIdPola /= pola
1887    else:       #'N'
1888        dIdPola = 0.0
1889    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1890    for iPO in dIdPO:
1891        dIdPO[iPO] /= POcorr
1892    dFdODF = {}
1893    dFdSA = [0,0,0]
1894    if pfx+'SHorder' in parmDict:
1895        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1896        for iSH in dFdODF:
1897            dFdODF[iSH] /= odfCor
1898        for i in range(3):
1899            dFdSA[i] /= odfCor
1900    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
1901       
1902def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1903    costh = cosd(refl[5]/2.)
1904    #crystallite size
1905    if calcControls[phfx+'SizeType'] == 'isotropic':
1906        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh)
1907    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1908        H = np.array(refl[:3])
1909        P = np.array(calcControls[phfx+'SizeAxis'])
1910        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1911        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']*costh)
1912        gam *= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)
1913    else:           #ellipsoidal crystallites - wrong not make sense
1914        H = np.array(refl[:3])
1915        gam += 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1916    #microstrain               
1917    if calcControls[phfx+'MustrainType'] == 'isotropic':
1918        gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
1919    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1920        H = np.array(refl[:3])
1921        P = np.array(calcControls[phfx+'MustrainAxis'])
1922        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1923        Si = parmDict[phfx+'Mustrain:0']
1924        Sa = parmDict[phfx+'Mustrain:1']
1925        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1926    else:       #generalized - P.W. Stephens model
1927        pwrs = calcControls[phfx+'MuPwrs']
1928        sum = 0
1929        for i,pwr in enumerate(pwrs):
1930            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1931        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
1932    return gam
1933       
1934def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1935    gamDict = {}
1936    costh = cosd(refl[5]/2.)
1937    tanth = tand(refl[5]/2.)
1938    #crystallite size derivatives
1939    if calcControls[phfx+'SizeType'] == 'isotropic':
1940        gamDict[phfx+'Size:0'] = -1.80*wave/(np.pi*costh)
1941    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1942        H = np.array(refl[:3])
1943        P = np.array(calcControls[phfx+'SizeAxis'])
1944        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1945        Si = parmDict[phfx+'Size:0']
1946        Sa = parmDict[phfx+'Size:1']
1947        gami = (1.80*wave/np.pi)/(Si*Sa)
1948        sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2)
1949        gam = gami*sqtrm/costh           
1950        gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)-gam/Si
1951        gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)-gam/Sa         
1952    else:           #ellipsoidal crystallites - do numerically? - not right not make sense
1953        H = np.array(refl[:3])
1954        gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1955    #microstrain derivatives               
1956    if calcControls[phfx+'MustrainType'] == 'isotropic':
1957        gamDict[phfx+'Mustrain:0'] =  0.018*tanth/np.pi           
1958    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1959        H = np.array(refl[:3])
1960        P = np.array(calcControls[phfx+'MustrainAxis'])
1961        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1962        Si = parmDict[phfx+'Mustrain:0']
1963        Sa = parmDict[phfx+'Mustrain:1']
1964        gami = 0.018*Si*Sa*tanth/np.pi
1965        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1966        gam = gami/sqtrm
1967        gamDict[phfx+'Mustrain:0'] = gam/Si-gami*Si*cosP**2/sqtrm**3
1968        gamDict[phfx+'Mustrain:1'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
1969    else:       #generalized - P.W. Stephens model
1970        pwrs = calcControls[phfx+'MuPwrs']
1971        const = 0.018*refl[4]**2*tanth
1972        for i,pwr in enumerate(pwrs):
1973            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1974    return gamDict
1975       
1976def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1977    h,k,l = refl[:3]
1978    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1979    d = np.sqrt(dsq)
1980    refl[4] = d
1981    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1982    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1983    if 'Bragg' in calcControls[hfx+'instType']:
1984        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1985            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1986    else:               #Debye-Scherrer - simple but maybe not right
1987        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1988    return pos
1989
1990def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1991    dpr = 180./np.pi
1992    h,k,l = refl[:3]
1993    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1994    dst = np.sqrt(dstsq)
1995    pos = refl[5]-parmDict[hfx+'Zero']
1996    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1997    dpdw = const*dst
1998    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1999    dpdA *= const*wave/(2.0*dst)
2000    dpdZ = 1.0
2001    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2002    if 'Bragg' in calcControls[hfx+'instType']:
2003        dpdSh = -4.*const*cosd(pos/2.0)
2004        dpdTr = -const*sind(pos)*100.0
2005        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2006    else:               #Debye-Scherrer - simple but maybe not right
2007        dpdXd = -const*cosd(pos)
2008        dpdYd = -const*sind(pos)
2009        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2010           
2011def GetHStrainShift(refl,SGData,phfx,parmDict):
2012    laue = SGData['SGLaue']
2013    uniq = SGData['SGUniq']
2014    h,k,l = refl[:3]
2015    if laue in ['m3','m3m']:
2016        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2017            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2018    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2019        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2020    elif laue in ['3R','3mR']:
2021        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2022    elif laue in ['4/m','4/mmm']:
2023        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2024    elif laue in ['mmm']:
2025        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2026    elif laue in ['2/m']:
2027        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2028        if uniq == 'a':
2029            Dij += parmDict[phfx+'D23']*k*l
2030        elif uniq == 'b':
2031            Dij += parmDict[phfx+'D13']*h*l
2032        elif uniq == 'c':
2033            Dij += parmDict[phfx+'D12']*h*k
2034    else:
2035        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2036            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2037    return Dij*refl[4]**2*tand(refl[5]/2.0)
2038           
2039def GetHStrainShiftDerv(refl,SGData,phfx):
2040    laue = SGData['SGLaue']
2041    uniq = SGData['SGUniq']
2042    h,k,l = refl[:3]
2043    if laue in ['m3','m3m']:
2044        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2045            phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2046    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2047        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2048    elif laue in ['3R','3mR']:
2049        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2050    elif laue in ['4/m','4/mmm']:
2051        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2052    elif laue in ['mmm']:
2053        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2054    elif laue in ['2/m']:
2055        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2056        if uniq == 'a':
2057            dDijDict[phfx+'D23'] = k*l
2058        elif uniq == 'b':
2059            dDijDict[phfx+'D13'] = h*l
2060        elif uniq == 'c':
2061            dDijDict[phfx+'D12'] = h*k
2062            names.append()
2063    else:
2064        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2065            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2066    for item in dDijDict:
2067        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
2068    return dDijDict
2069   
2070def GetFprime(controlDict,Histograms):
2071    FFtables = controlDict['FFtables']
2072    if not FFtables:
2073        return
2074    for histogram in Histograms:
2075        if 'PWDR' in histogram[:4]:
2076            Histogram = Histograms[histogram]
2077            hId = Histogram['hId']
2078            hfx = ':%d:'%(hId)
2079            keV = controlDict[hfx+'keV']
2080            for El in FFtables:
2081                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2082                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2083                FFtables[El][hfx+'FP'] = FP
2084                FFtables[El][hfx+'FPP'] = FPP               
2085           
2086def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2087   
2088    def GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse):
2089        U = parmDict[hfx+'U']
2090        V = parmDict[hfx+'V']
2091        W = parmDict[hfx+'W']
2092        X = parmDict[hfx+'X']
2093        Y = parmDict[hfx+'Y']
2094        tanPos = tand(refl[5]/2.0)
2095        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
2096        sig = max(0.001,sig)
2097        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
2098        gam = max(0.001,gam)
2099        return sig,gam
2100               
2101    hId = Histogram['hId']
2102    hfx = ':%d:'%(hId)
2103    bakType = calcControls[hfx+'bakType']
2104    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2105    yc = np.zeros_like(yb)
2106       
2107    if 'C' in calcControls[hfx+'histType']:   
2108        shl = max(parmDict[hfx+'SH/L'],0.002)
2109        Ka2 = False
2110        if hfx+'Lam1' in parmDict.keys():
2111            wave = parmDict[hfx+'Lam1']
2112            Ka2 = True
2113            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2114            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2115        else:
2116            wave = parmDict[hfx+'Lam']
2117    else:
2118        print 'TOF Undefined at present'
2119        raise ValueError
2120    for phase in Histogram['Reflection Lists']:
2121        refList = Histogram['Reflection Lists'][phase]
2122        Phase = Phases[phase]
2123        pId = Phase['pId']
2124        pfx = '%d::'%(pId)
2125        phfx = '%d:%d:'%(pId,hId)
2126        hfx = ':%d:'%(hId)
2127        SGData = Phase['General']['SGData']
2128        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2129        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2130        Vst = np.sqrt(nl.det(G))    #V*
2131        if 'Pawley' not in Phase['General']['Type']:
2132            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2133        sizeEllipse = []
2134        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
2135            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
2136        for refl in refList:
2137            if 'C' in calcControls[hfx+'histType']:
2138                h,k,l = refl[:3]
2139                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
2140                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
2141                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
2142                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
2143                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
2144                refl[13] *= Vst*Lorenz
2145                if 'Pawley' in Phase['General']['Type']:
2146                    try:
2147                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2148                    except KeyError:
2149#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2150                        continue
2151                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2152                iBeg = np.searchsorted(x,refl[5]-fmin)
2153                iFin = np.searchsorted(x,refl[5]+fmax)
2154                if not iBeg+iFin:       #peak below low limit - skip peak
2155                    continue
2156                elif not iBeg-iFin:     #peak above high limit - done
2157                    break
2158                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
2159                if Ka2:
2160                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2161                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2162                    iBeg = np.searchsorted(x,pos2-fmin)
2163                    iFin = np.searchsorted(x,pos2+fmax)
2164                    if not iBeg+iFin:       #peak below low limit - skip peak
2165                        continue
2166                    elif not iBeg-iFin:     #peak above high limit - done
2167                        return yc,yb
2168                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
2169            elif 'T' in calcControls[hfx+'histType']:
2170                print 'TOF Undefined at present'
2171                raise Exception    #no TOF yet
2172    return yc,yb
2173   
2174def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2175    for histogram in Histograms:
2176        if 'PWDR' in histogram[:4]:
2177            Histogram = Histograms[histogram]
2178            hId = Histogram['hId']
2179            hfx = ':%d:'%(hId)
2180            Limits = calcControls[hfx+'Limits']
2181            shl = max(parmDict[hfx+'SH/L'],0.002)
2182            Ka2 = False
2183            kRatio = 0.0
2184            if hfx+'Lam1' in parmDict.keys():
2185                Ka2 = True
2186                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2187                kRatio = parmDict[hfx+'I(L2)/I(L1)']
2188            x,y,w,yc,yb,yd = Histogram['Data']
2189            ymb = np.array(y-yb)
2190            ycmb = np.array(yc-yb)
2191            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
2192            xB = np.searchsorted(x,Limits[0])
2193            xF = np.searchsorted(x,Limits[1])
2194            refLists = Histogram['Reflection Lists']
2195            for phase in refLists:
2196                Phase = Phases[phase]
2197                pId = Phase['pId']
2198                phfx = '%d:%d:'%(pId,hId)
2199                refList = refLists[phase]
2200                sumFo = 0.0
2201                sumdF = 0.0
2202                sumFosq = 0.0
2203                sumdFsq = 0.0
2204                for refl in refList:
2205                    if 'C' in calcControls[hfx+'histType']:
2206                        yp = np.zeros_like(yb)
2207                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2208                        iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
2209                        iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
2210                        iFin2 = iFin
2211                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
2212                        if Ka2:
2213                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2214                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2215                            iBeg2 = np.searchsorted(x,pos2-fmin)
2216                            iFin2 = np.searchsorted(x,pos2+fmax)
2217                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
2218                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
2219                    elif 'T' in calcControls[hfx+'histType']:
2220                        print 'TOF Undefined at present'
2221                        raise Exception    #no TOF yet
2222                    Fo = np.sqrt(np.abs(refl[8]))
2223                    Fc = np.sqrt(np.abs(refl[9]))
2224                    sumFo += Fo
2225                    sumFosq += refl[8]**2
2226                    sumdF += np.abs(Fo-Fc)
2227                    sumdFsq += (refl[8]-refl[9])**2
2228                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2229                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2230                Histogram[phfx+'Nref'] = len(refList)
2231               
2232def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2233   
2234    def cellVaryDerv(pfx,SGData,dpdA): 
2235        if SGData['SGLaue'] in ['-1',]:
2236            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2237                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2238        elif SGData['SGLaue'] in ['2/m',]:
2239            if SGData['SGUniq'] == 'a':
2240                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2241            elif SGData['SGUniq'] == 'b':
2242                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2243            else:
2244                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2245        elif SGData['SGLaue'] in ['mmm',]:
2246            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2247        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2248#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2249            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2250        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2251#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2252            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2253        elif SGData['SGLaue'] in ['3R', '3mR']:
2254            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2255        elif SGData['SGLaue'] in ['m3m','m3']:
2256#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2257            return [[pfx+'A0',dpdA[0]]]
2258    # create a list of dependent variables and set up a dictionary to hold their derivatives
2259    dependentVars = G2mv.GetDependentVars()
2260    depDerivDict = {}
2261    for j in dependentVars:
2262        depDerivDict[j] = np.zeros(shape=(len(x)))
2263    #print 'dependent vars',dependentVars
2264    lenX = len(x)               
2265    hId = Histogram['hId']
2266    hfx = ':%d:'%(hId)
2267    bakType = calcControls[hfx+'bakType']
2268    dMdv = np.zeros(shape=(len(varylist),len(x)))
2269    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2270        dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2271        bBpos =varylist.index(hfx+'Back:0')
2272        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2273       
2274    if 'C' in calcControls[hfx+'histType']:   
2275        dx = x[1]-x[0]
2276        shl = max(parmDict[hfx+'SH/L'],0.002)
2277        Ka2 = False
2278        if hfx+'Lam1' in parmDict.keys():
2279            wave = parmDict[hfx+'Lam1']
2280            Ka2 = True
2281            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2282            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2283        else:
2284            wave = parmDict[hfx+'Lam']
2285    else:
2286        print 'TOF Undefined at present'
2287        raise ValueError
2288    for phase in Histogram['Reflection Lists']:
2289        refList = Histogram['Reflection Lists'][phase]
2290        Phase = Phases[phase]
2291        SGData = Phase['General']['SGData']
2292        pId = Phase['pId']
2293        pfx = '%d::'%(pId)
2294        phfx = '%d:%d:'%(pId,hId)
2295        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2296        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2297        if 'Pawley' not in Phase['General']['Type']:
2298            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2299        sizeEllipse = []
2300        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
2301            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
2302        for iref,refl in enumerate(refList):
2303            if 'C' in calcControls[hfx+'histType']:        #CW powder
2304                h,k,l = refl[:3]
2305                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2306                if 'Pawley' in Phase['General']['Type']:
2307                    try:
2308                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2309                    except KeyError:
2310#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2311                        continue
2312                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2313                iBeg = np.searchsorted(x,refl[5]-fmin)
2314                iFin = np.searchsorted(x,refl[5]+fmax)
2315                if not iBeg+iFin:       #peak below low limit - skip peak
2316                    continue
2317                elif not iBeg-iFin:     #peak above high limit - done
2318                    break
2319                pos = refl[5]
2320                tanth = tand(pos/2.0)
2321                costh = cosd(pos/2.0)
2322                dMdpk = np.zeros(shape=(6,len(x)))
2323                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2324                for i in range(1,5):
2325                    dMdpk[i][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
2326                dMdpk[0][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
2327                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
2328                if Ka2:
2329                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2330                    kdelt = int((pos2-refl[5])/dx)               
2331                    iBeg2 = min(lenX,iBeg+kdelt)
2332                    iFin2 = min(lenX,iFin+kdelt)
2333                    if iBeg2-iFin2:
2334                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2335                        for i in range(1,5):
2336                            dMdpk[i][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
2337                        dMdpk[0][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
2338                        dMdpk[5][iBeg2:iFin2] += 100.*dx*refl[13]*dMdipk2[0]
2339                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]}
2340                if 'Pawley' in Phase['General']['Type']:
2341                    try:
2342                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
2343                        dMdv[idx] = dervDict['int']/refl[9]
2344                        # Assuming Pawley variables not in constraints
2345                    except ValueError:
2346                        pass
2347                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2348                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2349                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2350                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2351                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2352                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2353                    hfx+'DisplaceY':[dpdY,'pos'],}
2354                for name in names:
2355                    item = names[name]
2356                    if name in varylist:
2357                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
2358                    elif name in dependentVars:
2359                        depDerivDict[name] += item[0]*dervDict[item[1]]
2360
2361                for iPO in dIdPO:
2362                    if iPO in varylist:
2363                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
2364                    elif iPO in dependentVars:
2365                        depDerivDict[iPO] = dIdPO[iPO]*dervDict['int']
2366
2367                for i,name in enumerate(['omega','chi','phi']):
2368                    aname = pfx+'SH '+name
2369                    if aname in varylist:
2370                        dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int']
2371                    elif aname in dependentVars:
2372                        depDerivDict[aname] += dFdSA[i]*dervDict['int']
2373                for iSH in dFdODF:
2374                    if iSH in varylist:
2375                        dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int']
2376                    elif iSH in dependentVars:
2377                        depDerivDict[iSH] += dFdODF[iSH]*dervDict['int']
2378                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2379                for name,dpdA in cellDervNames:
2380                    if name in varylist:
2381                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
2382                    elif name in dependentVars:
2383                        depDerivDict[name] += dpdA*dervDict['pos']
2384                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2385                for name in dDijDict:
2386                    if name in varylist:
2387                        dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos']
2388                    elif name in dependentVars:
2389                        depDerivDict[name] += dDijDict[name]*dervDict['pos']
2390                gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse)
2391                for name in gamDict:
2392                    if name in varylist:
2393                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']
2394                    elif name in dependentVars:
2395                        depDerivDict[name] += gamDict[name]*dervDict['gam']
2396                                               
2397            elif 'T' in calcControls[hfx+'histType']:
2398                print 'TOF Undefined at present'
2399                raise Exception    #no TOF yet
2400            #do atom derivatives -  for F,X & U so far             
2401            corr = dervDict['int']/refl[9]
2402            for name in varylist+dependentVars:
2403                try:
2404                    aname = name.split(pfx)[1][:2]
2405                    if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2406                except IndexError:
2407                    continue
2408                if name in varylist:
2409                    dMdv[varylist.index(name)] += dFdvDict[name][iref]*corr
2410                elif name in dependentVars:
2411                    depDerivDict[name] += dFdvDict[name][iref]*corr
2412    # now process derivatives in constraints
2413    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2414    return dMdv
2415
2416def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2417    parmdict.update(zip(varylist,values))
2418    G2mv.Dict2Map(parmdict,varylist)
2419    Histograms,Phases = HistoPhases
2420    dMdv = np.empty(0)
2421    for histogram in Histograms:
2422        if 'PWDR' in histogram[:4]:
2423            Histogram = Histograms[histogram]
2424            hId = Histogram['hId']
2425            hfx = ':%d:'%(hId)
2426            Limits = calcControls[hfx+'Limits']
2427            x,y,w,yc,yb,yd = Histogram['Data']
2428            xB = np.searchsorted(x,Limits[0])
2429            xF = np.searchsorted(x,Limits[1])
2430            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2431                varylist,Histogram,Phases,calcControls,pawleyLookup)
2432            if len(dMdv):
2433                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2434            else:
2435                dMdv = dMdvh
2436    return dMdv
2437
2438def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2439    parmdict.update(zip(varylist,values))
2440    Values2Dict(parmdict, varylist, values)
2441    G2mv.Dict2Map(parmdict,varylist)
2442    Histograms,Phases = HistoPhases
2443    M = np.empty(0)
2444    sumwYo = 0
2445    Nobs = 0
2446    for histogram in Histograms:
2447        if 'PWDR' in histogram[:4]:
2448            Histogram = Histograms[histogram]
2449            hId = Histogram['hId']
2450            hfx = ':%d:'%(hId)
2451            Limits = calcControls[hfx+'Limits']
2452            x,y,w,yc,yb,yd = Histogram['Data']
2453            yc *= 0.0                           #zero full calcd profiles
2454            yb *= 0.0
2455            yd *= 0.0
2456            xB = np.searchsorted(x,Limits[0])
2457            xF = np.searchsorted(x,Limits[1])
2458            Histogram['Nobs'] = xF-xB
2459            Nobs += Histogram['Nobs']
2460            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2461            sumwYo += Histogram['sumwYo']
2462            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2463                varylist,Histogram,Phases,calcControls,pawleyLookup)
2464            yc[xB:xF] += yb[xB:xF]
2465            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2466            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2467            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2468            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2469            M = np.concatenate((M,wdy))
2470    Histograms['sumwYo'] = sumwYo
2471    Histograms['Nobs'] = Nobs
2472    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2473    if dlg:
2474        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('wRp =',Rwp,'%'))[0]
2475        if not GoOn:
2476            parmDict['saved values'] = values
2477            raise Exception         #Abort!!
2478    return M
2479   
2480                   
2481def Refine(GPXfile,dlg):
2482    import cPickle
2483    import pytexture as ptx
2484    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2485   
2486    ShowBanner()
2487    varyList = []
2488    parmDict = {}
2489    calcControls = {}
2490    G2mv.InitVars()   
2491    Controls = GetControls(GPXfile)
2492    ShowControls(Controls)           
2493    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
2494    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2495    if not Phases:
2496        print ' *** ERROR - you have no histograms to refine! ***'
2497        print ' *** Refine aborted ***'
2498        raise Exception
2499    if not Histograms:
2500        print ' *** ERROR - you have no data to refine with! ***'
2501        print ' *** Refine aborted ***'
2502        raise Exception       
2503    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
2504    calcControls['Natoms'] = Natoms
2505    calcControls['FFtables'] = FFtables
2506    calcControls['BLtables'] = BLtables
2507    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2508    calcControls.update(controlDict)
2509    histVary,histDict,controlDict = GetHistogramData(Histograms)
2510    calcControls.update(controlDict)
2511    varyList = phaseVary+hapVary+histVary
2512    parmDict.update(phaseDict)
2513    parmDict.update(hapDict)
2514    parmDict.update(histDict)
2515    GetFprime(calcControls,Histograms)
2516    # do constraint processing
2517    try:
2518        groups,parmlist = G2mv.GroupConstraints(constrDict)
2519        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2520    except:
2521        print ' *** ERROR - your constraints are internally inconsistent ***'
2522        print ' *** Refine aborted ***'
2523        raise Exception
2524    # check to see which generated parameters are fully varied
2525    msg = G2mv.SetVaryFlags(varyList)
2526    if msg:
2527        print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
2528        print msg
2529        print ' *** Refine aborted ***'
2530        raise Exception       
2531    G2mv.Map2Dict(parmDict,varyList)
2532#    print G2mv.VarRemapShow(varyList)
2533
2534    while True:
2535        begin = time.time()
2536        values =  np.array(Dict2Values(parmDict, varyList))
2537        Ftol = Controls['min dM/M']
2538        Factor = Controls['shift factor']
2539        if Controls['deriv type'] == 'analytic':
2540            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2541                ftol=Ftol,col_deriv=True,factor=Factor,
2542                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2543            ncyc = int(result[2]['nfev']/2)               
2544        else:           #'numeric'
2545            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2546                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2547            ncyc = int(result[2]['nfev']/len(varyList))
2548#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
2549#        for item in table: print item,table[item]               #useful debug - are things shifting?
2550        runtime = time.time()-begin
2551        chisq = np.sum(result[2]['fvec']**2)
2552        Values2Dict(parmDict, varyList, result[0])
2553        G2mv.Dict2Map(parmDict,varyList)
2554        newCellDict = GetNewCellParms(parmDict,varyList)
2555        newAtomDict = ApplyXYZshifts(parmDict,varyList)
2556       
2557        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
2558        GOF = chisq/(Histograms['Nobs']-len(varyList))
2559        print '\n Refinement results:'
2560        print 135*'-'
2561        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
2562        print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2563        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
2564        try:
2565            covMatrix = result[1]*GOF
2566            sig = np.sqrt(np.diag(covMatrix))
2567            if np.any(np.isnan(sig)):
2568                print '*** Least squares aborted - some invalid esds possible ***'
2569#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
2570#            for item in table: print item,table[item]               #useful debug - are things shifting?
2571            break                   #refinement succeeded - finish up!
2572        except TypeError:          #result[1] is None on singular matrix
2573            print '**** Refinement failed - singular matrix ****'
2574            Ipvt = result[2]['ipvt']
2575            for i,ipvt in enumerate(Ipvt):
2576                if not np.sum(result[2]['fjac'],axis=1)[i]:
2577                    print 'Removing parameter: ',varyList[ipvt-1]
2578                    del(varyList[ipvt-1])
2579                    break
2580
2581#    print 'dependentParmList: ',G2mv.dependentParmList
2582#    print 'arrayList: ',G2mv.arrayList
2583#    print 'invarrayList: ',G2mv.invarrayList
2584#    print 'indParmList: ',G2mv.indParmList
2585#    print 'fixedDict: ',G2mv.fixedDict
2586#    print 'test1'
2587    GetFobsSq(Histograms,Phases,parmDict,calcControls)
2588#    print 'test2'
2589    sigDict = dict(zip(varyList,sig))
2590    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2591        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2592    # add the uncertainties into the esd dictionary (sigDict)
2593    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
2594    SetPhaseData(parmDict,sigDict,Phases,covData)
2595    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
2596    SetHistogramData(parmDict,sigDict,Histograms)
2597    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
2598    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
2599#for testing purposes!!!
2600#    file = open('structTestdata.dat','wb')
2601#    cPickle.dump(parmDict,file,1)
2602#    cPickle.dump(varyList,file,1)
2603#    for histogram in Histograms:
2604#        if 'PWDR' in histogram[:4]:
2605#            Histogram = Histograms[histogram]
2606#    cPickle.dump(Histogram,file,1)
2607#    cPickle.dump(Phases,file,1)
2608#    cPickle.dump(calcControls,file,1)
2609#    cPickle.dump(pawleyLookup,file,1)
2610#    file.close()
2611
2612def SeqRefine(GPXfile,dlg):
2613    import cPickle
2614    import pytexture as ptx
2615    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2616   
2617    ShowBanner()
2618    print ' Sequential Refinement'
2619    G2mv.InitVars()   
2620    Controls = GetControls(GPXfile)
2621    ShowControls(Controls)           
2622    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
2623    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2624    if not Phases:
2625        print ' *** ERROR - you have no histograms to refine! ***'
2626        print ' *** Refine aborted ***'
2627        raise Exception
2628    if not Histograms:
2629        print ' *** ERROR - you have no data to refine with! ***'
2630        print ' *** Refine aborted ***'
2631        raise Exception
2632    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
2633    if 'Seq Data' in Controls:
2634        histNames = Controls['Seq Data']
2635    else:
2636        histNames = GetHistogramNames(GPXfile,['PWDR',])
2637    if 'Reverse Seq' in Controls:
2638        if Controls['Reverse Seq']:
2639            histNames.reverse()
2640    SeqResult = {'histNames':histNames}
2641    makeBack = True
2642    for ihst,histogram in enumerate(histNames):
2643        ifPrint = False
2644        if dlg:
2645            dlg.SetTitle('Residual for histogram '+str(ihst))
2646        calcControls = {}
2647        calcControls['Natoms'] = Natoms
2648        calcControls['FFtables'] = FFtables
2649        calcControls['BLtables'] = BLtables
2650        varyList = []
2651        parmDict = {}
2652        Histo = {histogram:Histograms[histogram],}
2653        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
2654        calcControls.update(controlDict)
2655        histVary,histDict,controlDict = GetHistogramData(Histo,False)
2656        calcControls.update(controlDict)
2657        varyList = phaseVary+hapVary+histVary
2658        if not ihst:
2659            saveVaryList = varyList[:]
2660            for i,item in enumerate(saveVaryList):
2661                items = item.split(':')
2662                if items[1]:
2663                    items[1] = ''
2664                item = ':'.join(items)
2665                saveVaryList[i] = item
2666            SeqResult['varyList'] = saveVaryList
2667        else:
2668            newVaryList = varyList[:]
2669            for i,item in enumerate(newVaryList):
2670                items = item.split(':')
2671                if items[1]:
2672                    items[1] = ''
2673                item = ':'.join(items)
2674                newVaryList[i] = item
2675            if newVaryList != SeqResult['varyList']:
2676                print newVaryList
2677                print SeqResult['varyList']
2678                print '**** ERROR - variable list for this histogram does not match previous'
2679                raise Exception
2680        parmDict.update(phaseDict)
2681        parmDict.update(hapDict)
2682        parmDict.update(histDict)
2683        GetFprime(calcControls,Histo)
2684        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
2685        groups,parmlist = G2mv.GroupConstraints(constrDict)
2686        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2687        G2mv.Map2Dict(parmDict,varyList)
2688   
2689        while True:
2690            begin = time.time()
2691            values =  np.array(Dict2Values(parmDict, varyList))
2692            Ftol = Controls['min dM/M']
2693            Factor = Controls['shift factor']
2694            if Controls['deriv type'] == 'analytic':
2695                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2696                    ftol=Ftol,col_deriv=True,factor=Factor,
2697                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2698                ncyc = int(result[2]['nfev']/2)               
2699            else:           #'numeric'
2700                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2701                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2702                ncyc = int(result[2]['nfev']/len(varyList))
2703            runtime = time.time()-begin
2704            chisq = np.sum(result[2]['fvec']**2)
2705            Values2Dict(parmDict, varyList, result[0])
2706            G2mv.Dict2Map(parmDict,varyList)
2707            newCellDict = GetNewCellParms(parmDict,varyList)
2708            newAtomDict = ApplyXYZshifts(parmDict,varyList)
2709           
2710            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
2711            GOF = chisq/(Histo['Nobs']-len(varyList))
2712            print '\n Refinement results for histogram: v'+histogram
2713            print 135*'-'
2714            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
2715            print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2716            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
2717            try:
2718                covMatrix = result[1]*GOF
2719                sig = np.sqrt(np.diag(covMatrix))
2720                if np.any(np.isnan(sig)):
2721                    print '*** Least squares aborted - some invalid esds possible ***'
2722                    ifPrint = True
2723                break                   #refinement succeeded - finish up!
2724            except TypeError:          #result[1] is None on singular matrix
2725                print '**** Refinement failed - singular matrix ****'
2726                Ipvt = result[2]['ipvt']
2727                for i,ipvt in enumerate(Ipvt):
2728                    if not np.sum(result[2]['fjac'],axis=1)[i]:
2729                        print 'Removing parameter: ',varyList[ipvt-1]
2730                        del(varyList[ipvt-1])
2731                        break
2732   
2733        GetFobsSq(Histo,Phases,parmDict,calcControls)
2734        sigDict = dict(zip(varyList,sig))
2735        covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2736            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2737        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
2738        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
2739        SeqResult[histogram] = covData
2740        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
2741        makeBack = False
2742    SetSeqResult(GPXfile,Histograms,SeqResult)
2743
2744
2745def main():
2746    arg = sys.argv
2747    if len(arg) > 1:
2748        GPXfile = arg[1]
2749        if not ospath.exists(GPXfile):
2750            print 'ERROR - ',GPXfile," doesn't exist!"
2751            exit()
2752        GPXpath = ospath.dirname(arg[1])
2753        Refine(GPXfile,None)
2754    else:
2755        print 'ERROR - missing filename'
2756        exit()
2757    print "Done"
2758         
2759if __name__ == '__main__':
2760    main()
Note: See TracBrowser for help on using the repository browser.