source: trunk/GSASIIstruct.py @ 495

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

back to mine

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