source: trunk/GSASIIstruct.py @ 443

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

make sum powder profiles numpy arrays
change name UpdateBackgroundGrid? to UpdateBackground?
min size is 1 nanometer = 10A
implement diffuse scattering model for background
iso/aniso Size & Mustrain are now Size:i/a instead of 0/1

  • Property svn:keywords set to Date Author Revision URL Id
File size: 122.4 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2011-12-16 19:15:09 +0000 (Fri, 16 Dec 2011) $
4# $Author: vondreele $
5# $Revision: 443 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 443 2011-12-16 19:15:09Z 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+':i'] = hapData[item][1][0]
1071                    if hapData[item][2][0]:
1072                        hapVary.append(pfx+item+':i')
1073                    if hapData[item][0] == 'uniaxial':
1074                        controlDict[pfx+item+'Axis'] = hapData[item][3]
1075                        hapDict[pfx+item+':a'] = hapData[item][1][1]
1076                        if hapData[item][2][1]:
1077                            hapVary.append(pfx+item+':a')
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))+Zero
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+':i']
1260                    if item == 'Size':
1261                        hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1262                    if pfx+item+':i' in sigDict: 
1263                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
1264                    if hapData[item][0] == 'uniaxial':
1265                        hapData[item][1][1] = parmDict[pfx+item+':a']
1266                        if item == 'Size':
1267                            hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
1268                        if pfx+item+':a' in sigDict:
1269                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
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        Back = Background[0]
1290        Debye = Background[1]
1291        bakType,bakFlag = Back[:2]
1292        backVals = Back[3:]
1293        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1294        backDict = dict(zip(backNames,backVals))
1295        backVary = []
1296        if bakFlag:
1297            backVary = backNames
1298        backDict[':'+str(hId)+':nDebye'] = Debye['nDebye']
1299        debyeDict = {}
1300        debyeList = []
1301        for i in range(Debye['nDebye']):
1302            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
1303            debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1304            debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1305        debyeVary = []
1306        for item in debyeList:
1307            if item[1]:
1308                debyeVary.append(item[0])
1309        backDict.update(debyeDict)
1310        backVary += debyeVary   
1311        return bakType,backDict,backVary           
1312       
1313    def GetInstParms(hId,Inst):
1314        insVals,insFlags,insNames = Inst[1:4]
1315        dataType = insVals[0]
1316        instDict = {}
1317        insVary = []
1318        pfx = ':'+str(hId)+':'
1319        for i,flag in enumerate(insFlags):
1320            insName = pfx+insNames[i]
1321            instDict[insName] = insVals[i]
1322            if flag:
1323                insVary.append(insName)
1324        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
1325        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
1326        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1327        return dataType,instDict,insVary
1328       
1329    def GetSampleParms(hId,Sample):
1330        sampVary = []
1331        hfx = ':'+str(hId)+':'       
1332        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1333            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1334        Type = Sample['Type']
1335        if 'Bragg' in Type:             #Bragg-Brentano
1336            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
1337                sampDict[hfx+item] = Sample[item][0]
1338                if Sample[item][1]:
1339                    sampVary.append(hfx+item)
1340        elif 'Debye' in Type:        #Debye-Scherrer
1341            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1342                sampDict[hfx+item] = Sample[item][0]
1343                if Sample[item][1]:
1344                    sampVary.append(hfx+item)
1345        return Type,sampDict,sampVary
1346       
1347    def PrintBackground(Background):
1348        Back = Background[0]
1349        Debye = Background[1]
1350        print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
1351        line = ' Coefficients: '
1352        for i,back in enumerate(Back[3:]):
1353            line += '%10.3f'%(back)
1354            if i and not i%10:
1355                line += '\n'+15*' '
1356        print line
1357        if Debye['nDebye']:
1358            print '\n Debye diffuse scattering coefficients'
1359            parms = ['DebyeA','DebyeR','DebyeU']
1360            line = ' names :'
1361            for parm in parms:
1362                line += '%16s'%(parm)
1363            print line
1364            for j,term in enumerate(Debye['debyeTerms']):
1365                line = ' term'+'%2d'%(j)+':'
1366                for i in range(3):
1367                    line += '%10.4f %5s'%(term[2*i],bool(term[2*i+1]))                   
1368                print line
1369       
1370    def PrintInstParms(Inst):
1371        print '\n Instrument Parameters:'
1372        ptlbls = ' name  :'
1373        ptstr =  ' value :'
1374        varstr = ' refine:'
1375        instNames = Inst[3][1:]
1376        for i,name in enumerate(instNames):
1377            ptlbls += '%12s' % (name)
1378            ptstr += '%12.6f' % (Inst[1][i+1])
1379            if name in ['Lam1','Lam2','Azimuth']:
1380                varstr += 12*' '
1381            else:
1382                varstr += '%12s' % (str(bool(Inst[2][i+1])))
1383        print ptlbls
1384        print ptstr
1385        print varstr
1386       
1387    def PrintSampleParms(Sample):
1388        print '\n Sample Parameters:'
1389        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1390            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
1391        ptlbls = ' name  :'
1392        ptstr =  ' value :'
1393        varstr = ' refine:'
1394        if 'Bragg' in Sample['Type']:
1395            for item in ['Scale','Shift','Transparency']:
1396                ptlbls += '%14s'%(item)
1397                ptstr += '%14.4f'%(Sample[item][0])
1398                varstr += '%14s'%(str(bool(Sample[item][1])))
1399           
1400        elif 'Debye' in Type:        #Debye-Scherrer
1401            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1402                ptlbls += '%14s'%(item)
1403                ptstr += '%14.4f'%(Sample[item][0])
1404                varstr += '%14s'%(str(bool(Sample[item][1])))
1405
1406        print ptlbls
1407        print ptstr
1408        print varstr
1409       
1410
1411    histDict = {}
1412    histVary = []
1413    controlDict = {}
1414    for histogram in Histograms:
1415        Histogram = Histograms[histogram]
1416        hId = Histogram['hId']
1417        pfx = ':'+str(hId)+':'
1418        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1419       
1420        Background = Histogram['Background']
1421        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1422        controlDict[pfx+'bakType'] = Type
1423        histDict.update(bakDict)
1424        histVary += bakVary
1425       
1426        Inst = Histogram['Instrument Parameters']
1427        Type,instDict,insVary = GetInstParms(hId,Inst)
1428        controlDict[pfx+'histType'] = Type
1429        if pfx+'Lam1' in instDict:
1430            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1431        else:
1432            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
1433        histDict.update(instDict)
1434        histVary += insVary
1435       
1436        Sample = Histogram['Sample Parameters']
1437        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1438        controlDict[pfx+'instType'] = Type
1439        histDict.update(sampDict)
1440        histVary += sampVary
1441
1442        if Print: 
1443            print '\n Histogram: ',histogram,' histogram Id: ',hId
1444            print 135*'-'
1445            Units = {'C':' deg','T':' msec'}
1446            units = Units[controlDict[pfx+'histType'][2]]
1447            Limits = controlDict[pfx+'Limits']
1448            print ' Instrument type: ',Sample['Type']
1449            print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
1450            PrintSampleParms(Sample)
1451            PrintInstParms(Inst)
1452            PrintBackground(Background)
1453       
1454    return histVary,histDict,controlDict
1455   
1456def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
1457   
1458    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1459        Back = Background[0]
1460        Debye = Background[1]
1461        lenBack = len(Back[3:])
1462        backSig = [0 for i in range(lenBack+3*Debye['nDebye'])]
1463        for i in range(lenBack):
1464            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
1465            if pfx+'Back:'+str(i) in sigDict:
1466                backSig[i] = sigDict[pfx+'Back:'+str(i)]
1467        if Debye['nDebye']:
1468            for i in range(Debye['nDebye']):
1469                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
1470                for j,name in enumerate(names):
1471                    Debye['debyeTerms'][i][2*j] = parmDict[name]
1472                    if name in sigDict:
1473                        backSig[lenBack+3*i+j] = sigDict[name]           
1474        return backSig
1475       
1476    def SetInstParms(pfx,Inst,parmDict,sigDict):
1477        insVals,insFlags,insNames = Inst[1:4]
1478        instSig = [0 for i in range(len(insVals))]
1479        for i,flag in enumerate(insFlags):
1480            insName = pfx+insNames[i]
1481            insVals[i] = parmDict[insName]
1482            if insName in sigDict:
1483                instSig[i] = sigDict[insName]
1484        return instSig
1485       
1486    def SetSampleParms(pfx,Sample,parmDict,sigDict):
1487        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
1488            sampSig = [0 for i in range(3)]
1489            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
1490                Sample[item][0] = parmDict[pfx+item]
1491                if pfx+item in sigDict:
1492                    sampSig[i] = sigDict[pfx+item]
1493        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1494            sampSig = [0 for i in range(4)]
1495            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1496                Sample[item][0] = parmDict[pfx+item]
1497                if pfx+item in sigDict:
1498                    sampSig[i] = sigDict[pfx+item]
1499        return sampSig
1500       
1501    def PrintBackgroundSig(Background,backSig):
1502        Back = Background[0]
1503        Debye = Background[1]
1504        lenBack = len(Back[3:])
1505        print '\n Background function: ',Back[0]
1506        valstr = ' value : '
1507        sigstr = ' sig   : '
1508        for i,back in enumerate(Back[3:]):
1509            valstr += '%10.4f'%(back)
1510            if Back[1]:
1511                sigstr += '%10.4f'%(backSig[i])
1512            else:
1513                sigstr += 10*' '
1514        print valstr
1515        print sigstr
1516        if Debye['nDebye']:
1517            ifAny = False
1518            ptfmt = "%12.5f"
1519            names =  ' names :'
1520            ptstr =  ' values:'
1521            sigstr = ' esds  :'
1522            for item in sigDict:
1523                if 'Debye' in item:
1524                    ifAny = True
1525                    names += '%12s'%(item)
1526                    ptstr += ptfmt%(parmDict[item])
1527                    sigstr += ptfmt%(sigDict[item])
1528            if ifAny:
1529                print '\n Debye diffuse scattering coefficients'
1530                print names
1531                print ptstr
1532                print sigstr
1533       
1534    def PrintInstParmsSig(Inst,instSig):
1535        print '\n Instrument Parameters:'
1536        ptlbls = ' names :'
1537        ptstr =  ' value :'
1538        sigstr = ' sig   :'
1539        instNames = Inst[3][1:]
1540        for i,name in enumerate(instNames):
1541            ptlbls += '%12s' % (name)
1542            ptstr += '%12.6f' % (Inst[1][i+1])
1543            if instSig[i+1]:
1544                sigstr += '%12.6f' % (instSig[i+1])
1545            else:
1546                sigstr += 12*' '
1547        print ptlbls
1548        print ptstr
1549        print sigstr
1550       
1551    def PrintSampleParmsSig(Sample,sampleSig):
1552        print '\n Sample Parameters:'
1553        ptlbls = ' names :'
1554        ptstr =  ' values:'
1555        sigstr = ' sig   :'
1556        if 'Bragg' in Sample['Type']:
1557            for i,item in enumerate(['Scale','Shift','Transparency']):
1558                ptlbls += '%14s'%(item)
1559                ptstr += '%14.4f'%(Sample[item][0])
1560                if sampleSig[i]:
1561                    sigstr += '%14.4f'%(sampleSig[i])
1562                else:
1563                    sigstr += 14*' '
1564           
1565        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1566            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1567                ptlbls += '%14s'%(item)
1568                ptstr += '%14.4f'%(Sample[item][0])
1569                if sampleSig[i]:
1570                    sigstr += '%14.4f'%(sampleSig[i])
1571                else:
1572                    sigstr += 14*' '
1573
1574        print ptlbls
1575        print ptstr
1576        print sigstr
1577       
1578    for histogram in Histograms:
1579        if 'PWDR' in histogram:
1580            Histogram = Histograms[histogram]
1581            hId = Histogram['hId']
1582            pfx = ':'+str(hId)+':'
1583            Background = Histogram['Background']
1584            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1585           
1586            Inst = Histogram['Instrument Parameters']
1587            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1588       
1589            Sample = Histogram['Sample Parameters']
1590            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1591
1592            print '\n Histogram: ',histogram,' histogram Id: ',hId
1593            print 135*'-'
1594            print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
1595            if Print:
1596                print ' Instrument type: ',Sample['Type']
1597                PrintSampleParmsSig(Sample,sampSig)
1598                PrintInstParmsSig(Inst,instSig)
1599                PrintBackgroundSig(Background,backSig)
1600
1601def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict):
1602    Natoms = calcControls['Natoms'][pfx]
1603    Tdata = Natoms*[' ',]
1604    Mdata = np.zeros(Natoms)
1605    IAdata = Natoms*[' ',]
1606    Fdata = np.zeros(Natoms)
1607    FFdata = []
1608    BLdata = []
1609    Xdata = np.zeros((3,Natoms))
1610    dXdata = np.zeros((3,Natoms))
1611    Uisodata = np.zeros(Natoms)
1612    Uijdata = np.zeros((6,Natoms))
1613    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1614        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1615        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1616        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1617        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1618    for iatm in range(Natoms):
1619        for key in keys:
1620            parm = pfx+key+str(iatm)
1621            if parm in parmDict:
1622                keys[key][iatm] = parmDict[parm]
1623        FFdata.append(FFtables[Tdata[iatm]])
1624        BLdata.append(BLtables[Tdata[iatm]][1])
1625    return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1626   
1627def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1628    ''' Compute structure factors for all h,k,l for phase
1629    input:
1630        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1631        G:      reciprocal metric tensor
1632        pfx:    phase id string
1633        SGData: space group info. dictionary output from SpcGroup
1634        calcControls:
1635        ParmDict:
1636    puts result F^2 in each ref[8] in refList
1637    '''       
1638    twopi = 2.0*np.pi
1639    twopisq = 2.0*np.pi**2
1640    ast = np.sqrt(np.diag(G))
1641    Mast = twopisq*np.multiply.outer(ast,ast)
1642    FFtables = calcControls['FFtables']
1643    BLtables = calcControls['BLtables']
1644    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
1645    if 'N' in parmDict[hfx+'Type']:
1646        FP,FPP = G2el.BlenRes(BLdata,parmDict[hfx+'Lam'])
1647    else:
1648        FP = np.array([El[hfx+'FP'] for El in FFdata])
1649        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1650    maxPos = len(SGData['SGOps'])
1651    Uij = np.array(G2lat.U6toUij(Uijdata))
1652    bij = Mast*Uij.T
1653    for refl in refList:
1654        fbs = np.array([0,0])
1655        H = refl[:3]
1656        SQ = 1./(2.*refl[4])**2
1657        if 'N' in parmDict[hfx+'Type']:
1658            FF = np.array([El[1] for El in BLdata])
1659        else:       #'X'
1660            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1661        SQfactor = 4.0*SQ*twopisq
1662        Uniq = refl[11]
1663        phi = refl[12]
1664        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1665        sinp = np.sin(phase)
1666        cosp = np.cos(phase)
1667        occ = Mdata*Fdata/len(Uniq)
1668        biso = -SQfactor*Uisodata
1669        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1670        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1671        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1672        Tcorr = Tiso*Tuij
1673        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
1674        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
1675        if not SGData['SGInv']:
1676            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
1677            fbs = np.sum(np.sum(fb,axis=1),axis=1)
1678        fasq = fas**2
1679        fbsq = fbs**2        #imaginary
1680        refl[9] = np.sum(fasq)+np.sum(fbsq)
1681        refl[10] = atan2d(fbs[0],fas[0])
1682    return refList
1683   
1684def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1685    twopi = 2.0*np.pi
1686    twopisq = 2.0*np.pi**2
1687    ast = np.sqrt(np.diag(G))
1688    Mast = twopisq*np.multiply.outer(ast,ast)
1689    FFtables = calcControls['FFtables']
1690    BLtables = calcControls['BLtables']
1691    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
1692    if 'N' in parmDict[hfx+'Type']:
1693        FP = 0.
1694        FPP = 0.
1695    else:
1696        FP = np.array([El[hfx+'FP'] for El in FFdata])
1697        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1698    maxPos = len(SGData['SGOps'])       
1699    Uij = np.array(G2lat.U6toUij(Uijdata))
1700    bij = Mast*Uij.T
1701    dFdvDict = {}
1702    dFdfr = np.zeros((len(refList),len(Mdata)))
1703    dFdx = np.zeros((len(refList),len(Mdata),3))
1704    dFdui = np.zeros((len(refList),len(Mdata)))
1705    dFdua = np.zeros((len(refList),len(Mdata),6))
1706    for iref,refl in enumerate(refList):
1707        H = np.array(refl[:3])
1708        SQ = 1./(2.*refl[4])**2          # or (sin(theta)/lambda)**2
1709        if 'N' in parmDict[hfx+'Type']:
1710            FF = np.array([El[1] for El in BLdata])
1711        else:       #'X'
1712            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1713        SQfactor = 8.0*SQ*np.pi**2
1714        Uniq = refl[11]
1715        phi = refl[12]
1716        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
1717        sinp = np.sin(phase)
1718        cosp = np.cos(phase)
1719        occ = Mdata*Fdata/len(Uniq)
1720        biso = -SQfactor*Uisodata
1721        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1722#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1723        HbH = -np.inner(H,np.inner(bij,H))
1724        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1725        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1726        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1727        Tcorr = Tiso*Tuij
1728        fot = (FF+FP)*occ*Tcorr
1729        fotp = FPP*occ*Tcorr
1730        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1731        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1732       
1733        fas = np.sum(np.sum(fa,axis=1),axis=1)
1734        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1735        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1736        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1737        #sum below is over Uniq
1738        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
1739        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1740        dfadui = np.sum(-SQfactor*fa,axis=2)
1741        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1742        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
1743        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1744        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1745        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1746        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1747        if not SGData['SGInv']:
1748            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
1749            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
1750            dfbdui = np.sum(-SQfactor*fb,axis=2)
1751            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1752            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1753            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1754            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1755            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1756        #loop over atoms - each dict entry is list of derivatives for all the reflections
1757    for i in range(len(Mdata)):     
1758        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1759        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1760        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1761        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1762        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1763        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1764        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1765        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1766        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1767        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1768        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1769    return dFdvDict
1770       
1771def Dict2Values(parmdict, varylist):
1772    '''Use before call to leastsq to setup list of values for the parameters
1773    in parmdict, as selected by key in varylist'''
1774    return [parmdict[key] for key in varylist] 
1775   
1776def Values2Dict(parmdict, varylist, values):
1777    ''' Use after call to leastsq to update the parameter dictionary with
1778    values corresponding to keys in varylist'''
1779    parmdict.update(zip(varylist,values))
1780   
1781def GetNewCellParms(parmDict,varyList):
1782    newCellDict = {}
1783    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
1784    for item in varyList:
1785        keys = item.split(':')
1786        if keys[2] in Ddict:
1787            key = keys[0]+'::'+Ddict[keys[2]]
1788            parm = keys[0]+'::'+keys[2]
1789            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1790    return newCellDict
1791   
1792def ApplyXYZshifts(parmDict,varyList):
1793    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1794        input:
1795            parmDict - parameter dictionary
1796            varyList - list of variables
1797        returns:
1798            newAtomDict - dictitemionary of new atomic coordinate names & values;
1799                key is parameter shift name
1800    '''
1801    newAtomDict = {}
1802    for item in parmDict:
1803        if 'dA' in item:
1804            parm = ''.join(item.split('d'))
1805            parmDict[parm] += parmDict[item]
1806            newAtomDict[item] = [parm,parmDict[parm]]
1807    return newAtomDict
1808   
1809def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1810    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1811    odfCor = 1.0
1812    H = refl[:3]
1813    cell = G2lat.Gmat2cell(g)
1814    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1815    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1816    phi,beta = G2lat.CrsAng(H,cell,SGData)
1817    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1818    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1819    for item in SHnames:
1820        L,M,N = eval(item.strip('C'))
1821        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1822        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1823        Lnorm = G2lat.Lnorm(L)
1824        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1825    return odfCor
1826   
1827def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1828    FORPI = 12.5663706143592
1829    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1830    odfCor = 1.0
1831    dFdODF = {}
1832    dFdSA = [0,0,0]
1833    H = refl[:3]
1834    cell = G2lat.Gmat2cell(g)
1835    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1836    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1837    phi,beta = G2lat.CrsAng(H,cell,SGData)
1838    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
1839    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1840    for item in SHnames:
1841        L,M,N = eval(item.strip('C'))
1842        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1843        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1844        Lnorm = G2lat.Lnorm(L)
1845        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1846        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1847        for i in range(3):
1848            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1849    return odfCor,dFdODF,dFdSA
1850   
1851def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1852    odfCor = 1.0
1853    H = refl[:3]
1854    cell = G2lat.Gmat2cell(g)
1855    Sangl = [0.,0.,0.]
1856    if 'Bragg' in calcControls[hfx+'instType']:
1857        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1858        IFCoup = True
1859    else:
1860        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1861        IFCoup = False
1862    phi,beta = G2lat.CrsAng(H,cell,SGData)
1863    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1864    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1865    for item in SHnames:
1866        L,N = eval(item.strip('C'))
1867        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1868        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1869    return odfCor
1870   
1871def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1872    FORPI = 12.5663706143592
1873    odfCor = 1.0
1874    dFdODF = {}
1875    H = refl[:3]
1876    cell = G2lat.Gmat2cell(g)
1877    Sangl = [0.,0.,0.]
1878    if 'Bragg' in calcControls[hfx+'instType']:
1879        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1880        IFCoup = True
1881    else:
1882        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1883        IFCoup = False
1884    phi,beta = G2lat.CrsAng(H,cell,SGData)
1885    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1886    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1887    for item in SHnames:
1888        L,N = eval(item.strip('C'))
1889        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1890        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1891        dFdODF[phfx+item] = Kcsl*Lnorm
1892    return odfCor,dFdODF
1893   
1894def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1895    if calcControls[phfx+'poType'] == 'MD':
1896        MD = parmDict[phfx+'MD']
1897        MDAxis = calcControls[phfx+'MDAxis']
1898        sumMD = 0
1899        for H in refl[11]:           
1900            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1901            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1902            sumMD += A**3
1903        POcorr = sumMD/len(refl[11])
1904    else:   #spherical harmonics
1905        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1906    return POcorr
1907   
1908def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1909    POderv = {}
1910    if calcControls[phfx+'poType'] == 'MD':
1911        MD = parmDict[phfx+'MD']
1912        MDAxis = calcControls[phfx+'MDAxis']
1913        sumMD = 0
1914        sumdMD = 0
1915        for H in refl[11]:           
1916            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1917            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1918            sumMD += A**3
1919            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1920        POcorr = sumMD/len(refl[11])
1921        POderv[phfx+'MD'] = sumdMD/len(refl[11])
1922    else:   #spherical harmonics
1923        POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1924    return POcorr,POderv
1925   
1926def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1927    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1928    if 'X' in parmDict[hfx+'Type']:
1929        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1930    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1931    if pfx+'SHorder' in parmDict:
1932        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1933    refl[13] = Icorr       
1934   
1935def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1936    dIdsh = 1./parmDict[hfx+'Scale']
1937    dIdsp = 1./parmDict[phfx+'Scale']
1938    if 'X' in parmDict[hfx+'Type']:
1939        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1940        dIdPola /= pola
1941    else:       #'N'
1942        dIdPola = 0.0
1943    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1944    for iPO in dIdPO:
1945        dIdPO[iPO] /= POcorr
1946    dFdODF = {}
1947    dFdSA = [0,0,0]
1948    if pfx+'SHorder' in parmDict:
1949        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1950        for iSH in dFdODF:
1951            dFdODF[iSH] /= odfCor
1952        for i in range(3):
1953            dFdSA[i] /= odfCor
1954    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
1955       
1956def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict):
1957    costh = cosd(refl[5]/2.)
1958    #crystallite size
1959    if calcControls[phfx+'SizeType'] == 'isotropic':
1960        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
1961    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1962        H = np.array(refl[:3])
1963        P = np.array(calcControls[phfx+'SizeAxis'])
1964        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1965        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh)
1966        gam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2)
1967    else:           #ellipsoidal crystallites
1968        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1969        H = np.array(refl[:3])
1970        lenR = G2pwd.ellipseSize(H,Sij,GB)
1971        gam = 1.8*wave/(np.pi*costh*lenR)
1972    #microstrain               
1973    if calcControls[phfx+'MustrainType'] == 'isotropic':
1974        gam += 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
1975    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1976        H = np.array(refl[:3])
1977        P = np.array(calcControls[phfx+'MustrainAxis'])
1978        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1979        Si = parmDict[phfx+'Mustrain:i']
1980        Sa = parmDict[phfx+'Mustrain:a']
1981        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1982    else:       #generalized - P.W. Stephens model
1983        pwrs = calcControls[phfx+'MuPwrs']
1984        sum = 0
1985        for i,pwr in enumerate(pwrs):
1986            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1987        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
1988    return gam
1989       
1990def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
1991    gamDict = {}
1992    costh = cosd(refl[5]/2.)
1993    tanth = tand(refl[5]/2.)
1994    #crystallite size derivatives
1995    if calcControls[phfx+'SizeType'] == 'isotropic':
1996        gamDict[phfx+'Size:i'] = -1.80*wave/(np.pi*costh)
1997    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1998        H = np.array(refl[:3])
1999        P = np.array(calcControls[phfx+'SizeAxis'])
2000        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2001        Si = parmDict[phfx+'Size:i']
2002        Sa = parmDict[phfx+'Size:a']
2003        gami = (1.8*wave/np.pi)/(Si*Sa)
2004        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2005        gam = gami*sqtrm/costh           
2006        gamDict[phfx+'Size:i'] = gami*Si*cosP**2/(sqtrm*costh)-gam/Si
2007        gamDict[phfx+'Size:a'] = gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa         
2008    else:           #ellipsoidal crystallites
2009        const = 1.8*wave/(np.pi*costh)
2010        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2011        H = np.array(refl[:3])
2012        R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2013        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2014            gamDict[item] = -(const/R**2)*dRdS[i]
2015    #microstrain derivatives               
2016    if calcControls[phfx+'MustrainType'] == 'isotropic':
2017        gamDict[phfx+'Mustrain:i'] =  0.018*tanth/np.pi           
2018    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2019        H = np.array(refl[:3])
2020        P = np.array(calcControls[phfx+'MustrainAxis'])
2021        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2022        Si = parmDict[phfx+'Mustrain:i']
2023        Sa = parmDict[phfx+'Mustrain:a']
2024        gami = 0.018*Si*Sa*tanth/np.pi
2025        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2026        gam = gami/sqtrm
2027        gamDict[phfx+'Mustrain:i'] = gam/Si-gami*Si*cosP**2/sqtrm**3
2028        gamDict[phfx+'Mustrain:a'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
2029    else:       #generalized - P.W. Stephens model
2030        pwrs = calcControls[phfx+'MuPwrs']
2031        const = 0.018*refl[4]**2*tanth
2032        for i,pwr in enumerate(pwrs):
2033            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2034    return gamDict
2035       
2036def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
2037    h,k,l = refl[:3]
2038    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
2039    d = np.sqrt(dsq)
2040
2041    refl[4] = d
2042    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2043    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2044    if 'Bragg' in calcControls[hfx+'instType']:
2045        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2046            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2047    else:               #Debye-Scherrer - simple but maybe not right
2048        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2049    return pos
2050
2051def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
2052    dpr = 180./np.pi
2053    h,k,l = refl[:3]
2054    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2055    dst = np.sqrt(dstsq)
2056    pos = refl[5]-parmDict[hfx+'Zero']
2057    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2058    dpdw = const*dst
2059    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
2060    dpdA *= const*wave/(2.0*dst)
2061    dpdZ = 1.0
2062    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2063    if 'Bragg' in calcControls[hfx+'instType']:
2064        dpdSh = -4.*const*cosd(pos/2.0)
2065        dpdTr = -const*sind(pos)*100.0
2066        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2067    else:               #Debye-Scherrer - simple but maybe not right
2068        dpdXd = -const*cosd(pos)
2069        dpdYd = -const*sind(pos)
2070        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2071           
2072def GetHStrainShift(refl,SGData,phfx,parmDict):
2073    laue = SGData['SGLaue']
2074    uniq = SGData['SGUniq']
2075    h,k,l = refl[:3]
2076    if laue in ['m3','m3m']:
2077        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2078            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2079    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2080        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2081    elif laue in ['3R','3mR']:
2082        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2083    elif laue in ['4/m','4/mmm']:
2084        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2085    elif laue in ['mmm']:
2086        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2087    elif laue in ['2/m']:
2088        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2089        if uniq == 'a':
2090            Dij += parmDict[phfx+'D23']*k*l
2091        elif uniq == 'b':
2092            Dij += parmDict[phfx+'D13']*h*l
2093        elif uniq == 'c':
2094            Dij += parmDict[phfx+'D12']*h*k
2095    else:
2096        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2097            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2098    return Dij*refl[4]**2*tand(refl[5]/2.0)
2099           
2100def GetHStrainShiftDerv(refl,SGData,phfx):
2101    laue = SGData['SGLaue']
2102    uniq = SGData['SGUniq']
2103    h,k,l = refl[:3]
2104    if laue in ['m3','m3m']:
2105        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2106            phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2107    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2108        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2109    elif laue in ['3R','3mR']:
2110        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2111    elif laue in ['4/m','4/mmm']:
2112        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2113    elif laue in ['mmm']:
2114        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2115    elif laue in ['2/m']:
2116        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2117        if uniq == 'a':
2118            dDijDict[phfx+'D23'] = k*l
2119        elif uniq == 'b':
2120            dDijDict[phfx+'D13'] = h*l
2121        elif uniq == 'c':
2122            dDijDict[phfx+'D12'] = h*k
2123            names.append()
2124    else:
2125        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2126            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2127    for item in dDijDict:
2128        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
2129    return dDijDict
2130   
2131def GetFprime(controlDict,Histograms):
2132    FFtables = controlDict['FFtables']
2133    if not FFtables:
2134        return
2135    for histogram in Histograms:
2136        if 'PWDR' in histogram[:4]:
2137            Histogram = Histograms[histogram]
2138            hId = Histogram['hId']
2139            hfx = ':%d:'%(hId)
2140            keV = controlDict[hfx+'keV']
2141            for El in FFtables:
2142                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2143                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2144                FFtables[El][hfx+'FP'] = FP
2145                FFtables[El][hfx+'FPP'] = FPP               
2146           
2147def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2148   
2149    def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
2150        U = parmDict[hfx+'U']
2151        V = parmDict[hfx+'V']
2152        W = parmDict[hfx+'W']
2153        X = parmDict[hfx+'X']
2154        Y = parmDict[hfx+'Y']
2155        tanPos = tand(refl[5]/2.0)
2156        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
2157        sig = max(0.001,sig)
2158        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict) #save peak gamma
2159        gam = max(0.001,gam)
2160        return sig,gam
2161               
2162    hId = Histogram['hId']
2163    hfx = ':%d:'%(hId)
2164    bakType = calcControls[hfx+'bakType']
2165    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2166    yc = np.zeros_like(yb)
2167       
2168    if 'C' in calcControls[hfx+'histType']:   
2169        shl = max(parmDict[hfx+'SH/L'],0.002)
2170        Ka2 = False
2171        if hfx+'Lam1' in parmDict.keys():
2172            wave = parmDict[hfx+'Lam1']
2173            Ka2 = True
2174            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2175            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2176        else:
2177            wave = parmDict[hfx+'Lam']
2178    else:
2179        print 'TOF Undefined at present'
2180        raise ValueError
2181    for phase in Histogram['Reflection Lists']:
2182        refList = Histogram['Reflection Lists'][phase]
2183        Phase = Phases[phase]
2184        pId = Phase['pId']
2185        pfx = '%d::'%(pId)
2186        phfx = '%d:%d:'%(pId,hId)
2187        hfx = ':%d:'%(hId)
2188        SGData = Phase['General']['SGData']
2189        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2190        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2191        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2192        Vst = np.sqrt(nl.det(G))    #V*
2193        if 'Pawley' not in Phase['General']['Type']:
2194            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2195        for refl in refList:
2196            if 'C' in calcControls[hfx+'histType']:
2197                h,k,l = refl[:3]
2198                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
2199                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
2200                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
2201                refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
2202                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
2203                refl[13] *= Vst*Lorenz
2204                if 'Pawley' in Phase['General']['Type']:
2205                    try:
2206                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2207                    except KeyError:
2208#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2209                        continue
2210                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2211                iBeg = np.searchsorted(x,refl[5]-fmin)
2212                iFin = np.searchsorted(x,refl[5]+fmax)
2213                if not iBeg+iFin:       #peak below low limit - skip peak
2214                    continue
2215                elif not iBeg-iFin:     #peak above high limit - done
2216                    break
2217                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
2218                if Ka2:
2219                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2220                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2221                    iBeg = np.searchsorted(x,pos2-fmin)
2222                    iFin = np.searchsorted(x,pos2+fmax)
2223                    if not iBeg+iFin:       #peak below low limit - skip peak
2224                        continue
2225                    elif not iBeg-iFin:     #peak above high limit - done
2226                        return yc,yb
2227                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
2228            elif 'T' in calcControls[hfx+'histType']:
2229                print 'TOF Undefined at present'
2230                raise Exception    #no TOF yet
2231    return yc,yb
2232   
2233def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2234    for histogram in Histograms:
2235        if 'PWDR' in histogram[:4]:
2236            Histogram = Histograms[histogram]
2237            hId = Histogram['hId']
2238            hfx = ':%d:'%(hId)
2239            Limits = calcControls[hfx+'Limits']
2240            shl = max(parmDict[hfx+'SH/L'],0.002)
2241            Ka2 = False
2242            kRatio = 0.0
2243            if hfx+'Lam1' in parmDict.keys():
2244                Ka2 = True
2245                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2246                kRatio = parmDict[hfx+'I(L2)/I(L1)']
2247            x,y,w,yc,yb,yd = Histogram['Data']
2248            ymb = np.array(y-yb)
2249            ycmb = np.array(yc-yb)
2250            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
2251            xB = np.searchsorted(x,Limits[0])
2252            xF = np.searchsorted(x,Limits[1])
2253            refLists = Histogram['Reflection Lists']
2254            for phase in refLists:
2255                Phase = Phases[phase]
2256                pId = Phase['pId']
2257                phfx = '%d:%d:'%(pId,hId)
2258                refList = refLists[phase]
2259                sumFo = 0.0
2260                sumdF = 0.0
2261                sumFosq = 0.0
2262                sumdFsq = 0.0
2263                for refl in refList:
2264                    if 'C' in calcControls[hfx+'histType']:
2265                        yp = np.zeros_like(yb)
2266                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2267                        iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
2268                        iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
2269                        iFin2 = iFin
2270                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
2271                        if Ka2:
2272                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2273                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2274                            iBeg2 = np.searchsorted(x,pos2-fmin)
2275                            iFin2 = np.searchsorted(x,pos2+fmax)
2276                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
2277                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
2278                    elif 'T' in calcControls[hfx+'histType']:
2279                        print 'TOF Undefined at present'
2280                        raise Exception    #no TOF yet
2281                    Fo = np.sqrt(np.abs(refl[8]))
2282                    Fc = np.sqrt(np.abs(refl[9]))
2283                    sumFo += Fo
2284                    sumFosq += refl[8]**2
2285                    sumdF += np.abs(Fo-Fc)
2286                    sumdFsq += (refl[8]-refl[9])**2
2287                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2288                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2289                Histogram[phfx+'Nref'] = len(refList)
2290               
2291def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2292   
2293    def cellVaryDerv(pfx,SGData,dpdA): 
2294        if SGData['SGLaue'] in ['-1',]:
2295            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2296                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2297        elif SGData['SGLaue'] in ['2/m',]:
2298            if SGData['SGUniq'] == 'a':
2299                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2300            elif SGData['SGUniq'] == 'b':
2301                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2302            else:
2303                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2304        elif SGData['SGLaue'] in ['mmm',]:
2305            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2306        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2307#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2308            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2309        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2310#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2311            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2312        elif SGData['SGLaue'] in ['3R', '3mR']:
2313            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2314        elif SGData['SGLaue'] in ['m3m','m3']:
2315#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2316            return [[pfx+'A0',dpdA[0]]]
2317    # create a list of dependent variables and set up a dictionary to hold their derivatives
2318    dependentVars = G2mv.GetDependentVars()
2319    depDerivDict = {}
2320    for j in dependentVars:
2321        depDerivDict[j] = np.zeros(shape=(len(x)))
2322    #print 'dependent vars',dependentVars
2323    lenX = len(x)               
2324    hId = Histogram['hId']
2325    hfx = ':%d:'%(hId)
2326    bakType = calcControls[hfx+'bakType']
2327    dMdv = np.zeros(shape=(len(varylist),len(x)))
2328    dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2329    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2330        bBpos =varylist.index(hfx+'Back:0')
2331        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2332    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2333    for name in varylist:
2334        if 'Debye' in name:
2335            id = int(name.split(':')[-1])
2336            parm = name[:int(name.rindex(':'))]
2337            ip = names.index(parm)
2338            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2339    if 'C' in calcControls[hfx+'histType']:   
2340        dx = x[1]-x[0]
2341        shl = max(parmDict[hfx+'SH/L'],0.002)
2342        Ka2 = False
2343        if hfx+'Lam1' in parmDict.keys():
2344            wave = parmDict[hfx+'Lam1']
2345            Ka2 = True
2346            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2347            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2348        else:
2349            wave = parmDict[hfx+'Lam']
2350    else:
2351        print 'TOF Undefined at present'
2352        raise ValueError
2353    for phase in Histogram['Reflection Lists']:
2354        refList = Histogram['Reflection Lists'][phase]
2355        Phase = Phases[phase]
2356        SGData = Phase['General']['SGData']
2357        pId = Phase['pId']
2358        pfx = '%d::'%(pId)
2359        phfx = '%d:%d:'%(pId,hId)
2360        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2361        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2362        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2363        if 'Pawley' not in Phase['General']['Type']:
2364            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2365        for iref,refl in enumerate(refList):
2366            if 'C' in calcControls[hfx+'histType']:        #CW powder
2367                h,k,l = refl[:3]
2368                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2369                if 'Pawley' in Phase['General']['Type']:
2370                    try:
2371                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2372                    except KeyError:
2373#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2374                        continue
2375                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2376                iBeg = np.searchsorted(x,refl[5]-fmin)
2377                iFin = np.searchsorted(x,refl[5]+fmax)
2378                if not iBeg+iFin:       #peak below low limit - skip peak
2379                    continue
2380                elif not iBeg-iFin:     #peak above high limit - done
2381                    break
2382                pos = refl[5]
2383                tanth = tand(pos/2.0)
2384                costh = cosd(pos/2.0)
2385                dMdpk = np.zeros(shape=(6,len(x)))
2386                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2387                for i in range(1,5):
2388                    dMdpk[i][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
2389                dMdpk[0][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
2390                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
2391                if Ka2:
2392                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2393                    kdelt = int((pos2-refl[5])/dx)               
2394                    iBeg2 = min(lenX,iBeg+kdelt)
2395                    iFin2 = min(lenX,iFin+kdelt)
2396                    if iBeg2-iFin2:
2397                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2398                        for i in range(1,5):
2399                            dMdpk[i][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
2400                        dMdpk[0][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
2401                        dMdpk[5][iBeg2:iFin2] += 100.*dx*refl[13]*dMdipk2[0]
2402                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]}
2403                if 'Pawley' in Phase['General']['Type']:
2404                    try:
2405                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
2406                        dMdv[idx] = dervDict['int']/refl[9]
2407                        # Assuming Pawley variables not in constraints
2408                    except ValueError:
2409                        pass
2410                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2411                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2412                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2413                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2414                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2415                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2416                    hfx+'DisplaceY':[dpdY,'pos'],}
2417                for name in names:
2418                    item = names[name]
2419                    if name in varylist:
2420                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
2421                    elif name in dependentVars:
2422                        depDerivDict[name] += item[0]*dervDict[item[1]]
2423
2424                for iPO in dIdPO:
2425                    if iPO in varylist:
2426                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
2427                    elif iPO in dependentVars:
2428                        depDerivDict[iPO] = dIdPO[iPO]*dervDict['int']
2429
2430                for i,name in enumerate(['omega','chi','phi']):
2431                    aname = pfx+'SH '+name
2432                    if aname in varylist:
2433                        dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int']
2434                    elif aname in dependentVars:
2435                        depDerivDict[aname] += dFdSA[i]*dervDict['int']
2436                for iSH in dFdODF:
2437                    if iSH in varylist:
2438                        dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int']
2439                    elif iSH in dependentVars:
2440                        depDerivDict[iSH] += dFdODF[iSH]*dervDict['int']
2441                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2442                for name,dpdA in cellDervNames:
2443                    if name in varylist:
2444                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
2445                    elif name in dependentVars:
2446                        depDerivDict[name] += dpdA*dervDict['pos']
2447                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2448                for name in dDijDict:
2449                    if name in varylist:
2450                        dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos']
2451                    elif name in dependentVars:
2452                        depDerivDict[name] += dDijDict[name]*dervDict['pos']
2453                gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
2454                for name in gamDict:
2455                    if name in varylist:
2456                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']
2457                    elif name in dependentVars:
2458                        depDerivDict[name] += gamDict[name]*dervDict['gam']
2459                                               
2460            elif 'T' in calcControls[hfx+'histType']:
2461                print 'TOF Undefined at present'
2462                raise Exception    #no TOF yet
2463            #do atom derivatives -  for F,X & U so far             
2464            corr = dervDict['int']/refl[9]
2465            for name in varylist+dependentVars:
2466                try:
2467                    aname = name.split(pfx)[1][:2]
2468                    if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2469                except IndexError:
2470                    continue
2471                if name in varylist:
2472                    dMdv[varylist.index(name)] += dFdvDict[name][iref]*corr
2473                elif name in dependentVars:
2474                    depDerivDict[name] += dFdvDict[name][iref]*corr
2475    # now process derivatives in constraints
2476    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2477    return dMdv
2478
2479def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2480    parmdict.update(zip(varylist,values))
2481    G2mv.Dict2Map(parmdict,varylist)
2482    Histograms,Phases = HistoPhases
2483    nvar = len(varylist)
2484    dMdv = np.empty(0)
2485    for histogram in Histograms:
2486        if 'PWDR' in histogram[:4]:
2487            Histogram = Histograms[histogram]
2488            hId = Histogram['hId']
2489            hfx = ':%d:'%(hId)
2490            Limits = calcControls[hfx+'Limits']
2491            x,y,w,yc,yb,yd = Histogram['Data']
2492            xB = np.searchsorted(x,Limits[0])
2493            xF = np.searchsorted(x,Limits[1])
2494            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2495                varylist,Histogram,Phases,calcControls,pawleyLookup)
2496            if len(dMdv):
2497                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2498            else:
2499                dMdv = dMdvh
2500    return dMdv
2501
2502def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2503    parmdict.update(zip(varylist,values))
2504    Values2Dict(parmdict, varylist, values)
2505    G2mv.Dict2Map(parmdict,varylist)
2506    Histograms,Phases = HistoPhases
2507    M = np.empty(0)
2508    sumwYo = 0
2509    Nobs = 0
2510    for histogram in Histograms:
2511        if 'PWDR' in histogram[:4]:
2512            Histogram = Histograms[histogram]
2513            hId = Histogram['hId']
2514            hfx = ':%d:'%(hId)
2515            Limits = calcControls[hfx+'Limits']
2516            x,y,w,yc,yb,yd = Histogram['Data']
2517            yc *= 0.0                           #zero full calcd profiles
2518            yb *= 0.0
2519            yd *= 0.0
2520            xB = np.searchsorted(x,Limits[0])
2521            xF = np.searchsorted(x,Limits[1])
2522            Histogram['Nobs'] = xF-xB
2523            Nobs += Histogram['Nobs']
2524            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2525            sumwYo += Histogram['sumwYo']
2526            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2527                varylist,Histogram,Phases,calcControls,pawleyLookup)
2528            yc[xB:xF] += yb[xB:xF]
2529            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2530            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2531            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2532            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2533            M = np.concatenate((M,wdy))
2534    Histograms['sumwYo'] = sumwYo
2535    Histograms['Nobs'] = Nobs
2536    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2537    if dlg:
2538        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
2539        if not GoOn:
2540            parmDict['saved values'] = values
2541            raise Exception         #Abort!!
2542    return M
2543   
2544                   
2545def Refine(GPXfile,dlg):
2546    import cPickle
2547    import pytexture as ptx
2548    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2549   
2550    ShowBanner()
2551    varyList = []
2552    parmDict = {}
2553    calcControls = {}
2554    G2mv.InitVars()   
2555    Controls = GetControls(GPXfile)
2556    ShowControls(Controls)           
2557    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
2558    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2559    if not Phases:
2560        print ' *** ERROR - you have no histograms to refine! ***'
2561        print ' *** Refine aborted ***'
2562        raise Exception
2563    if not Histograms:
2564        print ' *** ERROR - you have no data to refine with! ***'
2565        print ' *** Refine aborted ***'
2566        raise Exception       
2567    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
2568    calcControls['Natoms'] = Natoms
2569    calcControls['FFtables'] = FFtables
2570    calcControls['BLtables'] = BLtables
2571    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2572    calcControls.update(controlDict)
2573    histVary,histDict,controlDict = GetHistogramData(Histograms)
2574    calcControls.update(controlDict)
2575    varyList = phaseVary+hapVary+histVary
2576    parmDict.update(phaseDict)
2577    parmDict.update(hapDict)
2578    parmDict.update(histDict)
2579    GetFprime(calcControls,Histograms)
2580    # do constraint processing
2581    try:
2582        groups,parmlist = G2mv.GroupConstraints(constrDict)
2583        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2584    except:
2585        print ' *** ERROR - your constraints are internally inconsistent ***'
2586        print ' *** Refine aborted ***'
2587        raise Exception
2588    # check to see which generated parameters are fully varied
2589    msg = G2mv.SetVaryFlags(varyList)
2590    if msg:
2591        print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
2592        print msg
2593        print ' *** Refine aborted ***'
2594        raise Exception       
2595    G2mv.Map2Dict(parmDict,varyList)
2596#    print G2mv.VarRemapShow(varyList)
2597
2598    while True:
2599        begin = time.time()
2600        values =  np.array(Dict2Values(parmDict, varyList))
2601        Ftol = Controls['min dM/M']
2602        Factor = Controls['shift factor']
2603        if Controls['deriv type'] == 'analytic':
2604            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2605                ftol=Ftol,col_deriv=True,factor=Factor,
2606                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2607            ncyc = int(result[2]['nfev']/2)               
2608        else:           #'numeric'
2609            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2610                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2611            ncyc = int(result[2]['nfev']/len(varyList))
2612#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
2613#        for item in table: print item,table[item]               #useful debug - are things shifting?
2614        runtime = time.time()-begin
2615        chisq = np.sum(result[2]['fvec']**2)
2616        Values2Dict(parmDict, varyList, result[0])
2617        G2mv.Dict2Map(parmDict,varyList)
2618       
2619        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
2620        GOF = chisq/(Histograms['Nobs']-len(varyList))
2621        print '\n Refinement results:'
2622        print 135*'-'
2623        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
2624        print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2625        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
2626        try:
2627            covMatrix = result[1]*GOF
2628            sig = np.sqrt(np.diag(covMatrix))
2629            if np.any(np.isnan(sig)):
2630                print '*** Least squares aborted - some invalid esds possible ***'
2631#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
2632#            for item in table: print item,table[item]               #useful debug - are things shifting?
2633            break                   #refinement succeeded - finish up!
2634        except TypeError:          #result[1] is None on singular matrix
2635            print '**** Refinement failed - singular matrix ****'
2636            Ipvt = result[2]['ipvt']
2637            for i,ipvt in enumerate(Ipvt):
2638                if not np.sum(result[2]['fjac'],axis=1)[i]:
2639                    print 'Removing parameter: ',varyList[ipvt-1]
2640                    del(varyList[ipvt-1])
2641                    break
2642
2643#    print 'dependentParmList: ',G2mv.dependentParmList
2644#    print 'arrayList: ',G2mv.arrayList
2645#    print 'invarrayList: ',G2mv.invarrayList
2646#    print 'indParmList: ',G2mv.indParmList
2647#    print 'fixedDict: ',G2mv.fixedDict
2648#    print 'test1'
2649    GetFobsSq(Histograms,Phases,parmDict,calcControls)
2650#    print 'test2'
2651    sigDict = dict(zip(varyList,sig))
2652    newCellDict = GetNewCellParms(parmDict,varyList)
2653    newAtomDict = ApplyXYZshifts(parmDict,varyList)
2654    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2655        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2656    # add the uncertainties into the esd dictionary (sigDict)
2657    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
2658    SetPhaseData(parmDict,sigDict,Phases,covData)
2659    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
2660    SetHistogramData(parmDict,sigDict,Histograms)
2661    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
2662    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
2663#for testing purposes!!!
2664#    file = open('structTestdata.dat','wb')
2665#    cPickle.dump(parmDict,file,1)
2666#    cPickle.dump(varyList,file,1)
2667#    for histogram in Histograms:
2668#        if 'PWDR' in histogram[:4]:
2669#            Histogram = Histograms[histogram]
2670#    cPickle.dump(Histogram,file,1)
2671#    cPickle.dump(Phases,file,1)
2672#    cPickle.dump(calcControls,file,1)
2673#    cPickle.dump(pawleyLookup,file,1)
2674#    file.close()
2675    return Rwp
2676
2677def SeqRefine(GPXfile,dlg):
2678    import cPickle
2679    import pytexture as ptx
2680    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2681   
2682    ShowBanner()
2683    print ' Sequential Refinement'
2684    G2mv.InitVars()   
2685    Controls = GetControls(GPXfile)
2686    ShowControls(Controls)           
2687    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
2688    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2689    if not Phases:
2690        print ' *** ERROR - you have no histograms to refine! ***'
2691        print ' *** Refine aborted ***'
2692        raise Exception
2693    if not Histograms:
2694        print ' *** ERROR - you have no data to refine with! ***'
2695        print ' *** Refine aborted ***'
2696        raise Exception
2697    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
2698    if 'Seq Data' in Controls:
2699        histNames = Controls['Seq Data']
2700    else:
2701        histNames = GetHistogramNames(GPXfile,['PWDR',])
2702    if 'Reverse Seq' in Controls:
2703        if Controls['Reverse Seq']:
2704            histNames.reverse()
2705    SeqResult = {'histNames':histNames}
2706    makeBack = True
2707    for ihst,histogram in enumerate(histNames):
2708        ifPrint = False
2709        if dlg:
2710            dlg.SetTitle('Residual for histogram '+str(ihst))
2711        calcControls = {}
2712        calcControls['Natoms'] = Natoms
2713        calcControls['FFtables'] = FFtables
2714        calcControls['BLtables'] = BLtables
2715        varyList = []
2716        parmDict = {}
2717        Histo = {histogram:Histograms[histogram],}
2718        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
2719        calcControls.update(controlDict)
2720        histVary,histDict,controlDict = GetHistogramData(Histo,False)
2721        calcControls.update(controlDict)
2722        varyList = phaseVary+hapVary+histVary
2723        if not ihst:
2724            saveVaryList = varyList[:]
2725            for i,item in enumerate(saveVaryList):
2726                items = item.split(':')
2727                if items[1]:
2728                    items[1] = ''
2729                item = ':'.join(items)
2730                saveVaryList[i] = item
2731            SeqResult['varyList'] = saveVaryList
2732        else:
2733            newVaryList = varyList[:]
2734            for i,item in enumerate(newVaryList):
2735                items = item.split(':')
2736                if items[1]:
2737                    items[1] = ''
2738                item = ':'.join(items)
2739                newVaryList[i] = item
2740            if newVaryList != SeqResult['varyList']:
2741                print newVaryList
2742                print SeqResult['varyList']
2743                print '**** ERROR - variable list for this histogram does not match previous'
2744                raise Exception
2745        parmDict.update(phaseDict)
2746        parmDict.update(hapDict)
2747        parmDict.update(histDict)
2748        GetFprime(calcControls,Histo)
2749        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
2750        groups,parmlist = G2mv.GroupConstraints(constrDict)
2751        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2752        G2mv.Map2Dict(parmDict,varyList)
2753   
2754        while True:
2755            begin = time.time()
2756            values =  np.array(Dict2Values(parmDict, varyList))
2757            Ftol = Controls['min dM/M']
2758            Factor = Controls['shift factor']
2759            if Controls['deriv type'] == 'analytic':
2760                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2761                    ftol=Ftol,col_deriv=True,factor=Factor,
2762                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2763                ncyc = int(result[2]['nfev']/2)               
2764            else:           #'numeric'
2765                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2766                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2767                ncyc = int(result[2]['nfev']/len(varyList))
2768            runtime = time.time()-begin
2769            chisq = np.sum(result[2]['fvec']**2)
2770            Values2Dict(parmDict, varyList, result[0])
2771            G2mv.Dict2Map(parmDict,varyList)
2772           
2773            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
2774            GOF = chisq/(Histo['Nobs']-len(varyList))
2775            print '\n Refinement results for histogram: v'+histogram
2776            print 135*'-'
2777            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
2778            print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2779            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
2780            try:
2781                covMatrix = result[1]*GOF
2782                sig = np.sqrt(np.diag(covMatrix))
2783                if np.any(np.isnan(sig)):
2784                    print '*** Least squares aborted - some invalid esds possible ***'
2785                    ifPrint = True
2786                break                   #refinement succeeded - finish up!
2787            except TypeError:          #result[1] is None on singular matrix
2788                print '**** Refinement failed - singular matrix ****'
2789                Ipvt = result[2]['ipvt']
2790                for i,ipvt in enumerate(Ipvt):
2791                    if not np.sum(result[2]['fjac'],axis=1)[i]:
2792                        print 'Removing parameter: ',varyList[ipvt-1]
2793                        del(varyList[ipvt-1])
2794                        break
2795   
2796        GetFobsSq(Histo,Phases,parmDict,calcControls)
2797        sigDict = dict(zip(varyList,sig))
2798        newCellDict = GetNewCellParms(parmDict,varyList)
2799        newAtomDict = ApplyXYZshifts(parmDict,varyList)
2800        covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2801            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2802        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
2803        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
2804        SeqResult[histogram] = covData
2805        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
2806        makeBack = False
2807    SetSeqResult(GPXfile,Histograms,SeqResult)
2808
2809
2810def main():
2811    arg = sys.argv
2812    if len(arg) > 1:
2813        GPXfile = arg[1]
2814        if not ospath.exists(GPXfile):
2815            print 'ERROR - ',GPXfile," doesn't exist!"
2816            exit()
2817        GPXpath = ospath.dirname(arg[1])
2818        Refine(GPXfile,None)
2819    else:
2820        print 'ERROR - missing filename'
2821        exit()
2822    print "Done"
2823         
2824if __name__ == '__main__':
2825    main()
Note: See TracBrowser for help on using the repository browser.