source: trunk/GSASIIstruct.py @ 453

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

refactor GSASIIgrid
new Hessian based least squares
new GSASIImath.py
work on focus issues

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