source: trunk/GSASIIstruct.py @ 367

Last change on this file since 367 was 367, checked in by vondreele, 12 years ago

make Factor lower limit 0.00001
  • Property svn:keywords set to Date Author Revision URL Id
File size: 67.3 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2011-09-08 20:27:13 +0000 (Thu, 08 Sep 2011) $
4# $Author: vondreele $
5# $Revision: 367 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 367 2011-09-08 20:27:13Z vondreele $
8########### SVN repository information ###################
9import sys
10import os
11import os.path as ospath
12import numpy as np
13import cPickle
14import time
15import math
16import wx
17import GSASIIpath
18import GSASIIElem as G2el
19import GSASIIlattice as G2lat
20import GSASIIspc as G2spc
21import GSASIIpwd as G2pwd
22import GSASIImapvars as G2mv
23import scipy.optimize as so
24
25sind = lambda x: np.sin(x*np.pi/180.)
26cosd = lambda x: np.cos(x*np.pi/180.)
27tand = lambda x: np.tan(x*np.pi/180.)
28asind = lambda x: 180.*np.arcsin(x)/np.pi
29
30
31def ShowBanner():
32    print 80*'*'
33    print '   General Structure Analysis System-II Crystal Structure Refinement'
34    print '     by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
35    print ' This product includes software developed by the UChicago Argonne, LLC,' 
36    print '            as Operator of Argonne National Laboratory.'
37    print 80*'*','\n'
38
39def GetControls(GPXfile):
40    ''' Returns dictionary of control items found in GSASII gpx file
41    input:
42        GPXfile = .gpx full file name
43    return:
44        Controls = dictionary of control items
45    '''
46    Controls = {'deriv type':'analytical','min dM/M':0.0001,'shift factor':1.}
47    file = open(GPXfile,'rb')
48    while True:
49        try:
50            data = cPickle.load(file)
51        except EOFError:
52            break
53        datum = data[0]
54        if datum[0] == 'Controls':
55            Controls.update(datum[1])
56    file.close()
57    return Controls
58   
59def ShowControls(Controls):
60    print ' Least squares controls:'
61    print ' Derivative type: ',Controls['deriv type']
62    print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
63    print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
64   
65def GetPhaseNames(GPXfile):
66    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
67    input:
68        GPXfile = gpx full file name
69    return:
70        PhaseNames = list of phase names
71    '''
72    file = open(GPXfile,'rb')
73    PhaseNames = []
74    while True:
75        try:
76            data = cPickle.load(file)
77        except EOFError:
78            break
79        datum = data[0]
80        if 'Phases' == datum[0]:
81            for datus in data[1:]:
82                PhaseNames.append(datus[0])
83    file.close()
84    return PhaseNames
85
86def GetAllPhaseData(GPXfile,PhaseName):
87    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
88    input:
89        GPXfile = gpx full file name
90        PhaseName = phase name
91    return:
92        phase dictionary
93    '''       
94    file = open(GPXfile,'rb')
95    General = {}
96    Atoms = []
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                if datus[0] == PhaseName:
106                    break
107    file.close()
108    return datus[1]
109   
110def GetHistogramNames(GPXfile):
111    ''' Returns a list of histogram names found in GSASII gpx file
112    input:
113        GPXfile = .gpx full file name
114    return:
115        HistogramNames = list of histogram names (types = PWDR & HKLF)
116    '''
117    file = open(GPXfile,'rb')
118    HistogramNames = []
119    while True:
120        try:
121            data = cPickle.load(file)
122        except EOFError:
123            break
124        datum = data[0]
125        if datum[0][:4] in ['PWDR','HKLF']:
126            HistogramNames.append(datum[0])
127    file.close()
128    return HistogramNames
129   
130def GetUsedHistogramsAndPhases(GPXfile):
131    ''' Returns all histograms that are found in any phase
132    and any phase that uses a histogram
133    input:
134        GPXfile = .gpx full file name
135    return:
136        Histograms = dictionary of histograms as {name:data,...}
137        Phases = dictionary of phases that use histograms
138    '''
139    phaseNames = GetPhaseNames(GPXfile)
140    phaseData = {}
141    for name in phaseNames: 
142        phaseData[name] =  GetAllPhaseData(GPXfile,name)
143    Histograms = {}
144    Phases = {}
145    pId = 0
146    hId = 0
147    for phase in phaseData:
148        Phase = phaseData[phase]
149        if Phase['Histograms']:
150            if phase not in Phases:
151                Phase['pId'] = pId
152                pId += 1
153                Phases[phase] = Phase
154            for hist in Phase['Histograms']:
155                if hist not in Histograms:
156                    if 'PWDR' in hist[:4]: 
157                        Histograms[hist] = GetPWDRdata(GPXfile,hist)
158                    elif 'HKLF' in hist[:4]:
159                        Histograms[hist] = GetHKLFdata(GPXfile,hist)
160                    #future restraint, etc. histograms here           
161                    Histograms[hist]['hId'] = hId
162                    hId += 1
163    return Histograms,Phases
164   
165def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases):
166    ''' Updates gpxfile from all histograms that are found in any phase
167    and any phase that used a histogram
168    input:
169        GPXfile = .gpx full file name
170        Histograms = dictionary of histograms as {name:data,...}
171        Phases = dictionary of phases that use histograms
172    '''
173                       
174    def GPXBackup(GPXfile):
175        import distutils.file_util as dfu
176        GPXpath,GPXname = ospath.split(GPXfile)
177        Name = ospath.splitext(GPXname)[0]
178        files = os.listdir(GPXpath)
179        last = 0
180        for name in files:
181            name = name.split('.')
182            if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
183                last = max(last,int(name[1].strip('bak'))+1)
184        GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
185        dfu.copy_file(GPXfile,GPXback)
186        return GPXback
187       
188    GPXback = GPXBackup(GPXfile)
189    print '\n',130*'-'
190    print 'Read from file:',GPXback
191    print 'Save to file  :',GPXfile
192    infile = open(GPXback,'rb')
193    outfile = open(GPXfile,'wb')
194    while True:
195        try:
196            data = cPickle.load(infile)
197        except EOFError:
198            break
199        datum = data[0]
200        print 'read: ',datum[0]
201        if datum[0] == 'Phases':
202            for iphase in range(len(data)):
203                if data[iphase][0] in Phases:
204                    phaseName = data[iphase][0]
205                    data[iphase][1] = Phases[phaseName]
206        try:
207            histogram = Histograms[datum[0]]
208            print 'found ',datum[0]
209            data[0][1][1] = histogram['Data']
210            for datus in data[1:]:
211                print '    read: ',datus[0]
212                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
213                    datus[1] = histogram[datus[0]]
214        except KeyError:
215            pass
216                               
217        cPickle.dump(data,outfile,1)
218    infile.close()
219    outfile.close()
220    print 'refinement save successful'
221                   
222def GetPWDRdata(GPXfile,PWDRname):
223    ''' Returns powder data from GSASII gpx file
224    input:
225        GPXfile = .gpx full file name
226        PWDRname = powder histogram name as obtained from GetHistogramNames
227    return:
228        PWDRdata = powder data dictionary with:
229            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
230       
231    '''
232    file = open(GPXfile,'rb')
233    PWDRdata = {}
234    while True:
235        try:
236            data = cPickle.load(file)
237        except EOFError:
238            break
239        datum = data[0]
240        if datum[0] == PWDRname:
241            PWDRdata['Data'] = datum[1][1]          #powder data arrays
242            PWDRdata[data[2][0]] = data[2][1]       #Limits
243            PWDRdata[data[3][0]] = data[3][1]       #Background
244            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
245            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
246            try:
247                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
248            except IndexError:
249                PWDRdata['Reflection lists'] = {}
250    file.close()
251    return PWDRdata
252   
253def GetHKLFdata(GPXfile,HKLFname):
254    ''' Returns single crystal data from GSASII gpx file
255    input:
256        GPXfile = .gpx full file name
257        HKLFname = single crystal histogram name as obtained from GetHistogramNames
258    return:
259        HKLFdata = single crystal data list of reflections: for each reflection:
260            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
261    '''
262    file = open(GPXfile,'rb')
263    HKLFdata = []
264    while True:
265        try:
266            data = cPickle.load(file)
267        except EOFError:
268            break
269        datum = data[0]
270        if datum[0] == HKLFname:
271            HKLFdata = datum[1:][0]
272    file.close()
273    return HKLFdata
274   
275def GetFFtable(General):
276    ''' returns a dictionary of form factor data for atom types found in General
277    input:
278        General = dictionary of phase info.; includes AtomTypes
279    return:
280        FFtable = dictionary of form factor data; key is atom type
281    '''
282    atomTypes = General['AtomTypes']
283    FFtable = {}
284    for El in atomTypes:
285        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
286        for item in FFs:
287            if item['Symbol'] == El:
288                FFtable[El] = item
289    return FFtable
290   
291def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
292    if SGLaue in ['-1','2/m','mmm']:
293        return                      #no Pawley symmetry required constraints
294    for i,varyI in enumerate(pawleyVary):
295        refI = int(varyI.split(':')[-1])
296        ih,ik,il = PawleyRef[refI][:3]
297        for varyJ in pawleyVary[0:i]:
298            refJ = int(varyJ.split(':')[-1])
299            jh,jk,jl = PawleyRef[refJ][:3]
300            if SGLaue in ['4/m','4/mmm']:
301                isum = ih**2+ik**2
302                jsum = jh**2+jk**2
303                if abs(il) == abs(jl) and isum == jsum:
304                    G2mv.StoreEquivalence(varyJ,(varyI,))
305            elif SGLaue in ['3R','3mR']:
306                isum = ih**2+ik**2+il**2
307                jsum = jh**2+jk**2*jl**2
308                isum2 = ih*ik+ih*il+ik*il
309                jsum2 = jh*jk+jh*jl+jk*jl
310                if isum == jsum and isum2 == jsum2:
311                    G2mv.StoreEquivalence(varyJ,(varyI,))
312            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
313                isum = ih**2+ik**2+ih*ik
314                jsum = jh**2+jk**2+jh*jk
315                if abs(il) == abs(jl) and isum == jsum:
316                    G2mv.StoreEquivalence(varyJ,(varyI,))
317            elif SGLaue in ['m3','m3m']:
318                isum = ih**2+ik**2+il**2
319                jsum = jh**2+jk**2+jl**2
320                if isum == jsum:
321                    G2mv.StoreEquivalence(varyJ,(varyI,))
322def cellVary(pfx,SGData): 
323    if SGData['SGLaue'] in ['-1',]:
324        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
325    elif SGData['SGLaue'] in ['2/m',]:
326        if SGData['SGUniq'] == 'a':
327            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
328        elif SGData['SGUniq'] == 'b':
329            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
330        else:
331            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
332    elif SGData['SGLaue'] in ['mmm',]:
333        return [pfx+'A0',pfx+'A1',pfx+'A2']
334    elif SGData['SGLaue'] in ['4/m','4/mmm']:
335        return [pfx+'A0',pfx+'A2']
336    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
337        return [pfx+'A0',pfx+'A2']
338    elif SGData['SGLaue'] in ['3R', '3mR']:
339        return [pfx+'A0',pfx+'A3']                       
340    elif SGData['SGLaue'] in ['m3m','m3']:
341        return [pfx+'A0']
342   
343       
344def GetPhaseData(PhaseData):
345   
346       
347    print ' Phases:'
348    phaseVary = []
349    phaseDict = {}
350    phaseConstr = {}
351    pawleyLookup = {}
352    for name in PhaseData:
353        General = PhaseData[name]['General']
354        pId = PhaseData[name]['pId']
355        pfx = str(pId)+'::'
356        Atoms = PhaseData[name]['Atoms']
357        try:
358            PawleyRef = PhaseData[name]['Pawley ref']
359        except KeyError:
360            PawleyRef = []
361        print '\n Phase name: ',General['Name']
362        SGData = General['SGData']
363        SGtext = G2spc.SGPrint(SGData)
364        for line in SGtext: print line
365        cell = General['Cell']
366        A = G2lat.cell2A(cell[1:7])
367        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]})
368        if cell[0]:
369            phaseVary = cellVary(pfx,SGData)
370        print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
371            ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
372            '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
373        if Atoms:
374            print '\n Atoms:'
375            line = '   name    type  refine?   x         y         z    '+ \
376                '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
377            if General['Type'] == 'magnetic':
378                line += '   Mx     My     Mz'
379            elif General['Type'] == 'macromolecular':
380                line = ' res no  residue  chain '+line
381            print line
382            if General['Type'] == 'nuclear':
383                print 135*'-'
384                for at in Atoms:
385                    line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
386                        '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
387                    if at[9] == 'I':
388                        line += '%8.4f'%(at[10])+48*' '
389                    else:
390                        line += 8*' '
391                        for i in range(6):
392                            line += '%8.4f'%(at[11+i])
393                    print line
394                    if 'X' in at[2]:
395                        xId,xCoef = G2spc.GetCSxinel(at[7])
396                    if 'U' in at[2]:
397                        uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
398#        elif General['Type'] == 'magnetic':
399#        elif General['Type'] == 'macromolecular':
400#       PWDR: harmonics texture, unit cell, atom parameters
401        elif PawleyRef:
402            pawleyVary = []
403            for i,refl in enumerate(PawleyRef):
404                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
405                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
406                if refl[5]:
407                    pawleyVary.append(pfx+'PWLref:'+str(i))
408            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
409            phaseVary += pawleyVary   
410               
411    return phaseVary,phaseDict,pawleyLookup
412   
413def SetPhaseData(parmDict,sigDict,Phases):
414   
415    def cellFill(pfx,SGData,parmDict,sigDict): 
416        if SGData['SGLaue'] in ['-1',]:
417            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
418                parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
419            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
420                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
421        elif SGData['SGLaue'] in ['2/m',]:
422            if SGData['SGUniq'] == 'a':
423                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
424                    parmDict[pfx+'A3'],0,0]
425                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
426                    sigDict[pfx+'A3'],0,0]
427            elif SGData['SGUniq'] == 'b':
428                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
429                    0,parmDict[pfx+'A4'],0]
430                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
431                    0,sigDict[pfx+'A4'],0]
432            else:
433                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
434                    0,0,parmDict[pfx+'A5']]
435                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
436                    0,0,sigDict[pfx+'A5']]
437        elif SGData['SGLaue'] in ['mmm',]:
438            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
439            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
440        elif SGData['SGLaue'] in ['4/m','4/mmm']:
441            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
442            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
443        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
444            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
445                parmDict[pfx+'A0'],0,0]
446            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
447        elif SGData['SGLaue'] in ['3R', '3mR']:
448            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
449                parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
450            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
451        elif SGData['SGLaue'] in ['m3m','m3']:
452            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
453            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
454        return A,sigA
455           
456    print '\n Phases:'
457    for phase in Phases:
458        print ' Result for phase: ',phase
459        Phase = Phases[phase]
460        General = Phase['General']
461        SGData = General['SGData']
462        cell = General['Cell']
463        pId = Phase['pId']
464        pfx = str(pId)+'::'
465        if cell[0]:
466            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
467            print ' Reciprocal metric tensor: '
468            ptfmt = "%15.9f"
469            names = ['A11','A22','A33','A12','A13','A23']
470            namstr = '  names :'
471            ptstr =  '  values:'
472            sigstr = '  esds  :'
473            for name,a,siga in zip(names,A,sigA):
474                namstr += '%15s'%(name)
475                ptstr += ptfmt%(a)
476                if siga:
477                    sigstr += ptfmt%(siga)
478                else:
479                    sigstr += 15*' '
480            print namstr
481            print ptstr
482            print sigstr
483            cell[1:7] = G2lat.A2cell(A)
484            cell[7] = G2lat.calc_V(A)
485            print ' New unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]), \
486                ' c = %.5f'%(cell[3]),' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]), \
487                ' gamma = %.3f'%(cell[6]),' volume = %.3f'%(cell[7])
488        if 'Pawley' in Phase['General']['Type']:
489            pawleyRef = Phase['Pawley ref']
490            for i,refl in enumerate(pawleyRef):
491                key = pfx+'PWLref:'+str(i)
492                refl[6] = abs(parmDict[key])        #suppress negative Fsq
493                if key in sigDict:
494                    refl[7] = sigDict[key]
495                else:
496                    refl[7] = 0
497
498def GetHistogramPhaseData(Phases,Histograms):
499   
500    def PrintSize(hapData):
501        line = '\n Size model    : '+hapData[0]
502        if hapData[0] in ['isotropic','uniaxial']:
503            line += ' equatorial:'+'%12.2f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
504            if hapData[0] == 'uniaxial':
505                line += ' axial:'+'%12.2f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
506            print line
507        else:
508            print line
509            Snames = ['S11','S22','S33','S12','S13','S23']
510            ptlbls = ' names :'
511            ptstr =  ' values:'
512            varstr = ' refine:'
513            for i,name in enumerate(Snames):
514                ptlbls += '%12s' % (name)
515                ptstr += '%12.6f' % (hapData[4][i])
516                varstr += '%12s' % (str(hapData[5][i]))
517            print ptlbls
518            print ptstr
519            print varstr
520       
521    def PrintMuStrain(hapData,SGData):
522        line = '\n Mustrain model: '+hapData[0]
523        if hapData[0] in ['isotropic','uniaxial']:
524            line += ' equatorial:'+'%12.4f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
525            if hapData[0] == 'uniaxial':
526                line += ' axial:'+'%12.4f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
527            print line
528        else:
529            print line
530            Snames = G2spc.MustrainNames(SGData)
531            ptlbls = ' names :'
532            ptstr =  ' values:'
533            varstr = ' refine:'
534            for i,name in enumerate(Snames):
535                ptlbls += '%12s' % (name)
536                ptstr += '%12.6f' % (hapData[4][i])
537                varstr += '%12s' % (str(hapData[5][i]))
538            print ptlbls
539            print ptstr
540            print varstr
541       
542   
543    hapDict = {}
544    hapVary = []
545    controlDict = {}
546    poType = {}
547    poAxes = {}
548    spAxes = {}
549    spType = {}
550   
551    for phase in Phases:
552        HistoPhase = Phases[phase]['Histograms']
553        SGData = Phases[phase]['General']['SGData']
554        cell = Phases[phase]['General']['Cell'][1:7]
555        A = G2lat.cell2A(cell)
556        pId = Phases[phase]['pId']
557        for histogram in HistoPhase:
558            print '\n Phase: ',phase,' in histogram: ',histogram
559            hapData = HistoPhase[histogram]
560            Histogram = Histograms[histogram]
561            hId = Histogram['hId']
562            limits = Histogram['Limits'][1]
563            inst = Histogram['Instrument Parameters']
564            inst = dict(zip(inst[3],inst[1]))
565            Zero = inst['Zero']
566            if 'C' in inst['Type']:
567                try:
568                    wave = inst['Lam']
569                except KeyError:
570                    wave = inst['Lam1']
571                dmin = wave/(2.0*sind(limits[1]/2.0))
572            pfx = str(pId)+':'+str(hId)+':'
573            for item in ['Scale','Extinction']:
574                hapDict[pfx+item] = hapData[item][0]
575                if hapData[item][1]:
576                    hapVary.append(pfx+item)
577            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
578            if hapData['Pref.Ori.'][0] == 'MD':
579                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
580                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
581                if hapData['Pref.Ori.'][2]:
582                    hapVary.append(pfx+'MD')
583            else:                           #'SH' spherical harmonics
584                for item in hapData['Pref.Ori.'][5]:
585                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
586                    if hapData['Pref.Ori.'][2]:
587                        hapVary.append(pfx+item)
588            for item in ['Mustrain','Size']:
589                controlDict[pfx+item+'Type'] = hapData[item][0]
590                if hapData[item][0] in ['isotropic','uniaxial']:
591                    hapDict[pfx+item+':0'] = hapData[item][1][0]
592                    if hapData[item][2][0]:
593                        hapVary.append(pfx+item+':0')
594                    if hapData[item][0] == 'uniaxial':
595                        controlDict[pfx+item+'Axis'] = hapData[item][3]
596                        hapDict[pfx+item+':1'] = hapData[item][1][1]
597                        if hapData[item][2][1]:
598                            hapVary.append(pfx+item+':1')
599                else:       #generalized for mustrain or ellipsoidal for size
600                    if item == 'Mustrain':
601                        names = G2spc.MustrainNames(SGData)
602                        pwrs = []
603                        for name in names:
604                            h,k,l = name[1:]
605                            pwrs.append([int(h),int(k),int(l)])
606                        controlDict[pfx+'MuPwrs'] = pwrs
607                    for i in range(len(hapData[item][4])):
608                        sfx = ':'+str(i)
609                        hapDict[pfx+item+sfx] = hapData[item][4][i]
610                        if hapData[item][5][i]:
611                            hapVary.append(pfx+item+sfx)
612                           
613            print '\n Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
614            print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
615            if hapData['Pref.Ori.'][0] == 'MD':
616                Ax = hapData['Pref.Ori.'][3]
617                print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
618                    ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])               
619            PrintSize(hapData['Size'])
620            PrintMuStrain(hapData['Mustrain'],SGData)
621            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
622            refList = []
623            for h,k,l,d in HKLd:
624                ext,mul,Uniq = G2spc.GenHKLf([h,k,l],SGData)
625                if ext:
626                    continue
627                if 'C' in inst['Type']:
628                    pos = 2.0*asind(wave/(2.0*d))
629                    if limits[0] < pos < limits[1]:
630                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,Uniq])
631                else:
632                    raise ValueError 
633            Histogram['Reflection Lists'][phase] = refList
634    return hapVary,hapDict,controlDict
635   
636def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms):
637   
638    def PrintSizeAndSig(hapData,sizeSig):
639        line = '\n Size model: '+hapData[0]
640        if hapData[0] in ['isotropic','uniaxial']:
641            line += ' equatorial:%12.2f'%(hapData[1][0])
642            if sizeSig[0][0]:
643                line += ', sig: %8.2f'%(sizeSig[0][0])
644            if hapData[0] == 'uniaxial':
645                line += ' axial:%12.2f'%(hapData[1][1])
646                if sizeSig[0][1]:
647                    line += ', sig: %8.2f'%(sizeSig[0][1])
648            print line
649        else:
650            print line
651            Snames = ['S11','S22','S33','S12','S13','S23']
652            ptlbls = ' name  :'
653            ptstr =  ' value :'
654            sigstr = ' sig   :'
655            for i,name in enumerate(Snames):
656                ptlbls += '%12s' % (name)
657                ptstr += '%12.6f' % (hapData[4][i])
658                if sizeSig[1][i]:
659                    sigstr += '%12.6f' % (sizeSig[1][i])
660                else:
661                    sigstr += 12*' '
662            print ptlbls
663            print ptstr
664            print sigstr
665       
666    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
667        line = '\n Mustrain model: '+hapData[0]
668        if hapData[0] in ['isotropic','uniaxial']:
669            line += ' equatorial:%12.4f'%(hapData[1][0])
670            if mustrainSig[0][0]:
671                line += ',sig: %12.4f'%(mustrainSig[0][0])
672            if hapData[0] == 'uniaxial':
673                line += ' axial:%12.4f'%(hapData[1][1])
674                if mustrainSig[0][1]:
675                     line += ',sig: %12.4f'%(mustrainSig[0][1])
676            print line
677        else:
678            print line
679            Snames = G2spc.MustrainNames(SGData)
680            ptlbls = ' name  :'
681            ptstr =  ' value :'
682            sigstr = ' sig   :'
683            for i,name in enumerate(Snames):
684                ptlbls += '%12s' % (name)
685                ptstr += '%12.6f' % (hapData[4][i])
686                if mustrainSig[1][i]:
687                    sigstr += '%12.6f' % (mustrainSig[1][i])
688                else:
689                    sigstr += 12*' '
690            print ptlbls
691            print ptstr
692            print sigstr
693       
694    for phase in Phases:
695        HistoPhase = Phases[phase]['Histograms']
696        SGData = Phases[phase]['General']['SGData']
697        pId = Phases[phase]['pId']
698        for histogram in HistoPhase:
699            print '\n Phase: ',phase,' in histogram: ',histogram
700            hapData = HistoPhase[histogram]
701            Histogram = Histograms[histogram]
702            hId = Histogram['hId']
703            pfx = str(pId)+':'+str(hId)+':'
704           
705            PhFrExtPOSig = {}
706            for item in ['Scale','Extinction']:
707                hapData[item][0] = parmDict[pfx+item]
708                if hapData[item][1]:
709                    PhFrExtPOSig[item] = sigDict[pfx+item]           
710            if hapData['Pref.Ori.'][0] == 'MD':
711                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
712                if hapData['Pref.Ori.'][2]:
713                    PhFrExtPOSig[item] = sigDict[pfx+item]
714            else:                           #'SH' spherical harmonics
715                for item in hapData['Pref.Ori.'][5]:
716                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
717                    if hapData['Pref.Ori.'][2]:
718                        PhFrExtPOSig[item] = sigDict[pfx+item]
719#            print '\n Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
720#            print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
721#            if hapData['Pref.Ori.'][0] == 'MD':
722#                Ax = hapData['Pref.Ori.'][3]
723#                print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
724#                    ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
725               
726            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
727                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]]}                 
728            for item in ['Mustrain','Size']:
729                if hapData[item][0] in ['isotropic','uniaxial']:                   
730                    hapData[item][1][0] = parmDict[pfx+item+':0']
731                    if item == 'Size':
732                        hapData[item][1][0] = min(10.,max(0.01,hapData[item][1][0]))
733                    if hapData[item][2][0]: 
734                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':0']
735                    if hapData[item][0] == 'uniaxial':
736                        hapData[item][1][1] = parmDict[pfx+item+':1']
737                        if item == 'Size':
738                            hapData[item][1][1] = min(10.,max(0.01,hapData[item][1][1]))                       
739                        if hapData[item][2][1]:
740                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':1']
741                else:       #generalized for mustrain or ellipsoidal for size
742                    for i in range(len(hapData[item][4])):
743                        sfx = ':'+str(i)
744                        hapData[item][4][i] = parmDict[pfx+item+sfx]
745                        if hapData[item][5][i]:
746                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
747                           
748            PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
749            PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
750   
751def GetHistogramData(Histograms):
752   
753    def GetBackgroundParms(hId,Background):
754        bakType,bakFlag = Background[:2]
755        backVals = Background[3:]
756        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
757        if bakFlag:                                 #returns backNames as varyList = backNames
758            return bakType,dict(zip(backNames,backVals)),backNames
759        else:                                       #no background varied; varyList = []
760            return bakType,dict(zip(backNames,backVals)),[]
761       
762    def GetInstParms(hId,Inst):
763        insVals,insFlags,insNames = Inst[1:4]
764        dataType = insVals[0]
765        instDict = {}
766        insVary = []
767        pfx = ':'+str(hId)+':'
768        for i,flag in enumerate(insFlags):
769            insName = pfx+insNames[i]
770            instDict[insName] = insVals[i]
771            if flag:
772                insVary.append(insName)
773        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.01)
774        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.01)
775        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
776        return dataType,instDict,insVary
777       
778    def GetSampleParms(hId,Sample):
779        sampVary = []
780        hfx = ':'+str(hId)+':'       
781        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius']}
782        Type = Sample['Type']
783        if 'Bragg' in Type:             #Bragg-Brentano
784            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
785                sampDict[hfx+item] = Sample[item][0]
786                if Sample[item][1]:
787                    sampVary.append(hfx+item)
788        elif 'Debye' in Type:        #Debye-Scherrer
789            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
790                sampDict[hfx+item] = Sample[item][0]
791                if Sample[item][1]:
792                    sampVary.append(hfx+item)
793        return Type,sampDict,sampVary
794       
795    def PrintBackground(Background):
796        print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
797        line = ' Coefficients: '
798        for back in Background[3:]:
799            line += '%10.3f'%(back)
800        print line
801       
802    def PrintInstParms(Inst):
803        print '\n Instrument Parameters:'
804        ptlbls = ' name  :'
805        ptstr =  ' value :'
806        varstr = ' refine:'
807        instNames = Inst[3][1:]
808        for i,name in enumerate(instNames):
809            ptlbls += '%12s' % (name)
810            ptstr += '%12.6f' % (Inst[1][i+1])
811            if name in ['Lam1','Lam2','Azimuth']:
812                varstr += 12*' '
813            else:
814                varstr += '%12s' % (str(bool(Inst[2][i+1])))
815        print ptlbls
816        print ptstr
817        print varstr
818       
819    def PrintSampleParms(Sample):
820        print '\n Sample Parameters:'
821        ptlbls = ' name  :'
822        ptstr =  ' value :'
823        varstr = ' refine:'
824        if 'Bragg' in Sample['Type']:
825            for item in ['Scale','Shift','Transparency']:
826                ptlbls += '%14s'%(item)
827                ptstr += '%14.4f'%(Sample[item][0])
828                varstr += '%14s'%(str(bool(Sample[item][1])))
829           
830        elif 'Debye' in Type:        #Debye-Scherrer
831            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
832                ptlbls += '%14s'%(item)
833                ptstr += '%14.4f'%(Sample[item][0])
834                varstr += '%14s'%(str(bool(Sample[item][1])))
835
836        print ptlbls
837        print ptstr
838        print varstr
839       
840
841    histDict = {}
842    histVary = []
843    controlDict = {}
844    for histogram in Histograms:
845        Histogram = Histograms[histogram]
846        hId = Histogram['hId']
847        pfx = ':'+str(hId)+':'
848        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
849       
850        Background = Histogram['Background'][0]
851        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
852        controlDict[pfx+'bakType'] = Type
853        histDict.update(bakDict)
854        histVary += bakVary
855       
856        Inst = Histogram['Instrument Parameters']
857        Type,instDict,insVary = GetInstParms(hId,Inst)
858        controlDict[pfx+'histType'] = Type
859        histDict.update(instDict)
860        histVary += insVary
861       
862        Sample = Histogram['Sample Parameters']
863        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
864        controlDict[pfx+'instType'] = Type
865        histDict.update(sampDict)
866        histVary += sampVary
867
868        print '\n Histogram: ',histogram,' histogram Id: ',hId
869        print 135*'-'
870        Units = {'C':' deg','T':' msec'}
871        units = Units[controlDict[pfx+'histType'][2]]
872        Limits = controlDict[pfx+'Limits']
873        print ' Instrument type: ',Sample['Type']
874        print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
875        PrintSampleParms(Sample)
876        PrintInstParms(Inst)
877        PrintBackground(Background)
878       
879    return histVary,histDict,controlDict
880   
881def SetHistogramData(parmDict,sigDict,Histograms):
882   
883    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
884        lenBack = len(Background[3:])
885        backSig = [0 for i in range(lenBack)]
886        for i in range(lenBack):
887            Background[3+i] = parmDict[pfx+'Back:'+str(i)]
888            if Background[1]:
889                backSig[i] = sigDict[pfx+'Back:'+str(i)]
890        return backSig
891       
892    def SetInstParms(pfx,Inst,parmDict,sigDict):
893        insVals,insFlags,insNames = Inst[1:4]
894        instSig = [0 for i in range(len(insVals))]
895        for i,flag in enumerate(insFlags):
896            insName = pfx+insNames[i]
897            insVals[i] = parmDict[insName]
898            if flag:
899                instSig[i] = sigDict[insName]
900        return instSig
901       
902    def SetSampleParms(pfx,Sample,parmDict,sigDict):
903        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
904            sampSig = [0 for i in range(3)]
905            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
906                Sample[item][0] = parmDict[pfx+item]
907                if Sample[item][1]:
908                    sampSig[i] = sigDict[pfx+item]
909        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
910            sampSig = [0 for i in range(4)]
911            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
912                Sample[item][0] = parmDict[pfx+item]
913                if Sample[item][1]:
914                    sampSig[i] = sigDict[pfx+item]
915        return sampSig
916       
917    def PrintBackgroundSig(Background,backSig):
918        print '\n Background function: ',Background[0]
919        valstr = ' value : '
920        sigstr = ' sig   : '
921        for i,back in enumerate(Background[3:]):
922            valstr += '%10.4f'%(back)
923            if Background[1]:
924                sigstr += '%10.4f'%(backSig[i])
925            else:
926                sigstr += 10*' '
927        print valstr
928        print sigstr
929       
930    def PrintInstParmsSig(Inst,instSig):
931        print '\n Instrument Parameters:'
932        ptlbls = ' names :'
933        ptstr =  ' value :'
934        sigstr = ' sig   :'
935        instNames = Inst[3][1:]
936        for i,name in enumerate(instNames):
937            ptlbls += '%12s' % (name)
938            ptstr += '%12.6f' % (Inst[1][i+1])
939            if instSig[i+1]:
940                sigstr += '%12.6f' % (instSig[i+1])
941            else:
942                sigstr += 12*' '
943        print ptlbls
944        print ptstr
945        print sigstr
946       
947    def PrintSampleParmsSig(Sample,sampleSig):
948        print '\n Sample Parameters:'
949        ptlbls = ' names :'
950        ptstr =  ' values:'
951        sigstr = ' sig   :'
952        if 'Bragg' in Sample['Type']:
953            for i,item in enumerate(['Scale','Shift','Transparency']):
954                ptlbls += '%14s'%(item)
955                ptstr += '%14.4f'%(Sample[item][0])
956                sigstr += '%14.4f'%(sampleSig[i])
957           
958        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
959            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
960                ptlbls += '%14s'%(item)
961                ptstr += '%14.4f'%(Sample[item][0])
962                sigstr += '%14.4f'%(sampleSig[i])
963
964        print ptlbls
965        print ptstr
966        print sigstr
967       
968    for histogram in Histograms:
969        if 'PWDR' in histogram:
970            Histogram = Histograms[histogram]
971            hId = Histogram['hId']
972            pfx = ':'+str(hId)+':'
973            Background = Histogram['Background'][0]
974            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
975           
976            Inst = Histogram['Instrument Parameters']
977            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
978       
979            Sample = Histogram['Sample Parameters']
980            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
981
982            print '\n Histogram: ',histogram,' histogram Id: ',hId
983            print 135*'-'
984            print ' Instrument type: ',Sample['Type']
985            PrintSampleParmsSig(Sample,sampSig)
986            PrintInstParmsSig(Inst,instSig)
987            PrintBackgroundSig(Background,backSig)
988
989
990def SetupSFcalc(General,Atoms):
991    ''' setup data for use in StructureFactor; mostly rearranging arrays
992    input:
993        General = dictionary of general phase info.
994        Atoms = list of atom parameters
995    returns:
996        G = reciprocal metric tensor
997        StrData = list of arrays; one entry per atom:
998            T = atom types
999            R = refinement flags, e.g. 'FXU'
1000            F = site fractions
1001            X = atom coordinates as numpy array
1002            M = atom multiplicities
1003            IA = isotropic/anisotropic thermal flags
1004            Uiso = isotropic thermal parameters if IA = 'I'; else = 0
1005            Uij = numpy array of 6 anisotropic thermal parameters if IA='A'; else = zeros
1006    '''           
1007    SGData = General['SGData']
1008    Cell = General['Cell']
1009    G,g = G2lat.cell2Gmat(Cell[1:7])        #skip refine & volume; get recip & real metric tensors
1010    cx,ct,cs,cia = General['AtomPtrs']
1011    X = [];F = [];T = [];IA = [];Uiso = [];Uij = [];R = [];M = []
1012    for atom in Atoms:
1013        T.append(atom[ct])
1014        R.append(atom[ct+1])
1015        F.append(atom[cx+3])
1016        X.append(np.array(atom[cx:cx+3]))
1017        M.append(atom[cia-1])
1018        IA.append(atom[cia])
1019        Uiso.append(atom[cia+1])
1020        Uij.append(np.array(atom[cia+2:cia+8]))
1021    return G,[T,R,np.array(F),np.array(X),np.array(M),IA,np.array(Uiso),np.array(Uij)]
1022           
1023def StructureFactor(H,G,SGData,StrData,FFtable):
1024    ''' Compute structure factor for a single h,k,l
1025    H = np.array(h,k,l)
1026    G = reciprocal metric tensor
1027    SGData = space group info. dictionary output from SpcGroup
1028    StrData = [
1029        [atom types],
1030        [refinement flags],
1031        [site fractions],
1032        np.array(coordinates),
1033        [multiplicities],
1034        [I/A flag],
1035        [Uiso],
1036        np.array(Uij)]
1037    FFtable = dictionary of form factor coeff. for atom types used in StrData
1038    '''
1039    twopi = 2.0*math.pi
1040    twopisq = 2.0*math.pi**2
1041    SQ = G2lat.calc_rDsq2(H,G)/4.0          # SQ = (sin(theta)/lambda)**2
1042    SQfactor = 8.0*SQ*math.pi**2
1043    print 'SQ',SQfactor
1044    FF = {}
1045    for type in FFtable.keys():
1046        FF[type] = G2el.ScatFac(FFtable[type],SQ)
1047    print 'FF',FF
1048    iabsnt,mulp,Uniq,Phs = G2spc.GenHKL(H,SGData,False)
1049    fa = [0,0,0]        #real
1050    fb = [0,0,0]        #imaginary
1051    if not iabsnt:
1052        phase = twopi*np.inner(Uniq,StrData[3])       #2pi[hx+ky+lz] for each atom in each equiv. hkl
1053        sinp = np.sin(phase)
1054        cosp = np.cos(phase)
1055        occ = StrData[2]*StrData[4]/mulp
1056        Tiso = np.multiply(StrData[6],SQfactor)
1057        Tuij = np.multiply(-SQfactor,np.inner(H,np.inner(G2spc.Uij2U(StrData[7]),H)))
1058        print 'sinp',sinp
1059        print 'cosp',cosp
1060        print 'occ',occ
1061        print 'Tiso',Tiso
1062        print 'Tuij',Tuij
1063    else:
1064        print 'Absent'
1065       
1066def Dict2Values(parmdict, varylist):
1067    '''Use before call to leastsq to setup list of values for the parameters
1068    in parmdict, as selected by key in varylist'''
1069    return [parmdict[key] for key in varylist] 
1070   
1071def Values2Dict(parmdict, varylist, values):
1072    ''' Use after call to leastsq to update the parameter dictionary with
1073    values corresponding to keys in varylist'''
1074    parmdict.update(zip(varylist,values))
1075   
1076def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1077   
1078    def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1079        costh = cosd(refl[5]/2.)
1080        #crystallite size
1081        if calcControls[phfx+'SizeType'] == 'isotropic':
1082            gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh)
1083        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1084            H = np.array(refl[:3])
1085            P = np.array(calcControls[phfx+'SizeAxis'])
1086            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1087            gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']*costh)
1088            gam *= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)
1089        else:           #ellipsoidal crystallites - wrong not make sense
1090            H = np.array(refl[:3])
1091            gam += 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1092        #microstrain               
1093        if calcControls[phfx+'MustrainType'] == 'isotropic':
1094            gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
1095        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1096            H = np.array(refl[:3])
1097            P = np.array(calcControls[phfx+'MustrainAxis'])
1098            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1099            Si = parmDict[phfx+'Mustrain:0']
1100            Sa = parmDict[phfx+'Mustrain:1']
1101            gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1102        else:       #generalized - P.W. Stephens model
1103            pwrs = calcControls[phfx+'MuPwrs']
1104            sum = 0
1105            for i,pwr in enumerate(pwrs):
1106                sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1107            gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
1108        return gam
1109       
1110    def GetIntensityCorr(refl,phfx,hfx,calcControls,parmDict):
1111        Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1112        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1113       
1114        return Icorr
1115       
1116    def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1117        h,k,l = refl[:3]
1118        dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1119        d = np.sqrt(dsq)
1120        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1121        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1122        if 'Bragg' in calcControls[hfx+'instType']:
1123            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1124                1.e-7*parmDict[hfx+'Transparency']*sind(pos))            #trans(=1/mueff) in Angstroms
1125        else:               #Debye-Scherrer - simple but maybe not right
1126            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1127        return pos
1128   
1129    def GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse):
1130        U = parmDict[hfx+'U']
1131        V = parmDict[hfx+'V']
1132        W = parmDict[hfx+'W']
1133        X = parmDict[hfx+'X']
1134        Y = parmDict[hfx+'Y']
1135        tanPos = tand(refl[5]/2.0)
1136        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
1137        sig = max(0.001,sig)
1138        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
1139        gam = max(0.01,gam)
1140        return sig,gam
1141               
1142    hId = Histogram['hId']
1143    hfx = ':%d:'%(hId)
1144    bakType = calcControls[hfx+'bakType']
1145    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1146    yc = np.zeros_like(yb)
1147       
1148    if 'C' in calcControls[hfx+'histType']:   
1149        shl = max(parmDict[hfx+'SH/L'],0.0005)
1150        Ka2 = False
1151        if hfx+'Lam1' in parmDict.keys():
1152            wave = parmDict[hfx+'Lam1']
1153            Ka2 = True
1154            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1155            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1156        else:
1157            wave = parmDict[hfx+'Lam']
1158    else:
1159        print 'TOF Undefined at present'
1160        raise ValueError
1161    for phase in Histogram['Reflection Lists']:
1162        refList = Histogram['Reflection Lists'][phase]
1163        Phase = Phases[phase]
1164        pId = Phase['pId']
1165        pfx = '%d::'%(pId)
1166        phfx = '%d:%d:'%(pId,hId)
1167        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1168        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1169        sizeEllipse = []
1170        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1171            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1172        for refl in refList:
1173            if 'C' in calcControls[hfx+'histType']:
1174                h,k,l = refl[:3]
1175                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1176                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
1177                Icorr = GetIntensityCorr(refl,phfx,hfx,calcControls,parmDict)
1178                if 'Pawley' in Phase['General']['Type']:
1179                    try:
1180                        refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
1181                    except KeyError:
1182                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1183                        continue
1184                else:
1185                    raise ValueError       #wants strctrfacr calc here
1186                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1187                iBeg = np.searchsorted(x,refl[5]-fmin)
1188                iFin = np.searchsorted(x,refl[5]+fmax)
1189                if not iBeg+iFin:       #peak below low limit - skip peak
1190                    continue
1191                elif not iBeg-iFin:     #peak above high limit - done
1192                    return yc,yb
1193                yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1194                if Ka2:
1195                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1196                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
1197                    iBeg = np.searchsorted(x,pos2-fmin)
1198                    iFin = np.searchsorted(x,pos2+fmax)
1199                    if not iBeg+iFin:       #peak below low limit - skip peak
1200                        continue
1201                    elif not iBeg-iFin:     #peak above high limit - done
1202                        return yc,yb
1203                    yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
1204            else:
1205                raise ValueError
1206    return yc,yb   
1207           
1208def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1209   
1210    def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1211        gamDict = {}
1212        costh = cosd(refl[5]/2.)
1213        tanth = tand(refl[5]/2.)
1214        #crystallite size derivatives
1215        if calcControls[phfx+'SizeType'] == 'isotropic':
1216            gamDict[phfx+'Size:0'] = 1.80*wave/(np.pi*costh)
1217        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1218            H = np.array(refl[:3])
1219            P = np.array(calcControls[phfx+'SizeAxis'])
1220            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1221            Si = parmDict[phfx+'Size:0']
1222            Sa = parmDict[phfx+'Size:1']
1223            gami = (1.80*wave/np.pi)/(Si*Sa)
1224            sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2)
1225            gam = gami*sqtrm/costh           
1226            gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)-gam/Si
1227            gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)-gam/Sa         
1228        else:           #ellipsoidal crystallites - do numerically? - not right not make sense
1229            H = np.array(refl[:3])
1230            gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1231        #microstrain derivatives               
1232        if calcControls[phfx+'MustrainType'] == 'isotropic':
1233            gamDict[phfx+'Mustrain:0'] =  0.018*tanth/np.pi           
1234        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1235            H = np.array(refl[:3])
1236            P = np.array(calcControls[phfx+'MustrainAxis'])
1237            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1238            Si = parmDict[phfx+'Mustrain:0']
1239            Sa = parmDict[phfx+'Mustrain:1']
1240            gami = 0.018*Si*Sa*tanth/np.pi
1241            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1242            gam = gami/sqtrm
1243            gamDict[phfx+'Mustrain:0'] = gam/Si-gami*Si*cosP**2/sqtrm**3
1244            gamDict[phfx+'Mustrain:1'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
1245        else:       #generalized - P.W. Stephens model
1246            pwrs = calcControls[phfx+'MuPwrs']
1247            const = 0.018*refl[4]**2*tanth
1248            for i,pwr in enumerate(pwrs):
1249                gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1250        return gamDict
1251       
1252    def GetIntensityDerv(refl,phfx,hfx,calcControls,parmDict):
1253        Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1254        pola,dpdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1255        dIdpola = Icorr*dpdPola
1256        Icorr *= pola
1257        dIdsh = Icorr/parmDict[hfx+'Scale']
1258        dIdsp = Icorr/parmDict[phfx+'Scale']
1259       
1260        return Icorr,dIdsh,dIdsp,dIdpola
1261       
1262    def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1263        dpr = 180./np.pi
1264        h,k,l = refl[:3]
1265        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1266        dst = np.sqrt(dstsq)
1267        pos = refl[5]
1268        const = dpr/np.sqrt(1.0-wave*dst/4.0)
1269        dpdw = const*dst
1270        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1271        dpdA *= const*wave/(2.0*dst)
1272        dpdZ = 1.0
1273        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1274        if 'Bragg' in calcControls[hfx+'instType']:
1275            dpdSh = -4.*const*cosd(pos/2.0)
1276            dpdTr = -1.e-7*const*sind(pos)
1277            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1278        else:               #Debye-Scherrer - simple but maybe not right
1279            dpdXd = -const*cosd(pos)
1280            dpdYd = -const*sind(pos)
1281            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1282           
1283    def cellVaryDerv(pfx,SGData,dpdA): 
1284        if SGData['SGLaue'] in ['-1',]:
1285            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1286                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1287        elif SGData['SGLaue'] in ['2/m',]:
1288            if SGData['SGUniq'] == 'a':
1289                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1290            elif SGData['SGUniq'] == 'b':
1291                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1292            else:
1293                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1294        elif SGData['SGLaue'] in ['mmm',]:
1295            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1296        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1297            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
1298        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1299            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
1300        elif SGData['SGLaue'] in ['3R', '3mR']:
1301            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1302        elif SGData['SGLaue'] in ['m3m','m3']:
1303            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
1304   
1305    lenX = len(x)               
1306    hId = Histogram['hId']
1307    hfx = ':%d:'%(hId)
1308    bakType = calcControls[hfx+'bakType']
1309    dMdv = np.zeros(shape=(len(varylist),len(x)))
1310    if hfx+'Back:0' in varylist:
1311        dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1312        bBpos =varylist.index(hfx+'Back:0')
1313        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1314       
1315    if 'C' in calcControls[hfx+'histType']:   
1316        dx = x[1]-x[0]
1317        shl = max(parmDict[hfx+'SH/L'],0.002)
1318        Ka2 = False
1319        if hfx+'Lam1' in parmDict.keys():
1320            wave = parmDict[hfx+'Lam1']
1321            Ka2 = True
1322            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1323            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1324        else:
1325            wave = parmDict[hfx+'Lam']
1326    else:
1327        print 'TOF Undefined at present'
1328        raise ValueError
1329    for phase in Histogram['Reflection Lists']:
1330        refList = Histogram['Reflection Lists'][phase]
1331        Phase = Phases[phase]
1332        SGData = Phase['General']['SGData']
1333        pId = Phase['pId']
1334        pfx = '%d::'%(pId)
1335        phfx = '%d:%d:'%(pId,hId)
1336        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1337        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1338        sizeEllipse = []
1339        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1340            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1341        for iref,refl in enumerate(refList):
1342            if 'C' in calcControls[hfx+'histType']:
1343                Icorr,dIdsh,dIdsp,dIdpola = GetIntensityDerv(refl,phfx,hfx,calcControls,parmDict)
1344                hkl = refl[:3]
1345                pos = refl[5]
1346                tanth = tand(pos/2.0)
1347                costh = cosd(pos/2.0)
1348                dsdU = tanth**2
1349                dsdV = tanth
1350                dsdW = 1.0
1351                dgdX = 1.0/costh
1352                dgdY = tanth
1353                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1354                iBeg = np.searchsorted(x,refl[5]-fmin)
1355                iFin = np.searchsorted(x,refl[5]+fmax)
1356                dMdpk = np.zeros(shape=(6,len(x)))
1357                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
1358                for i in range(1,5):
1359                    dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i]
1360                dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*dMdipk[0]
1361                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1362                if Ka2:
1363                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1364                    kdelt = int((pos2-refl[5])/dx)               
1365                    iBeg = min(lenX,iBeg+kdelt)
1366                    iFin = min(lenX,iFin+kdelt)
1367                    if iBeg-iFin:
1368                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])
1369                        for i in range(1,5):
1370                            dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]
1371                        dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
1372                        dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
1373                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*Icorr*refl[8]}
1374                try:
1375                    idx = varylist.index(pfx+'PWLref:'+str(iref))
1376                    dMdv[idx] = dervDict['int']
1377                except ValueError:
1378                    pass
1379                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1380                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1381                    hfx+'U':[dsdU,'sig'],hfx+'V':[dsdV,'sig'],hfx+'W':[dsdW,'sig'],
1382                    hfx+'X':[dgdX,'gam'],hfx+'Y':[dgdY,'gam'],hfx+'SH/L':[1.0,'shl'],
1383                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1384                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1385                    hfx+'DisplaceY':[dpdY,'pos'],}
1386                for name in names:
1387                    if name in varylist:
1388                        item = names[name]
1389                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
1390                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1391                for name,dpdA in cellDervNames:
1392                    if name in varylist:
1393                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
1394                gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse)
1395                for name in gamDict:
1396                    if name in varylist:
1397                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']                       
1398            else:
1399                raise ValueError
1400           
1401    return dMdv   
1402                   
1403def Refine(GPXfile):
1404   
1405    def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
1406        parmdict.update(zip(varylist,values))
1407        G2mv.Dict2Map(parmDict)
1408        Histograms,Phases = HistoPhases
1409        dMdv = np.empty(0)
1410        for histogram in Histograms:
1411            if 'PWDR' in histogram[:4]:
1412                Histogram = Histograms[histogram]
1413                hId = Histogram['hId']
1414                hfx = ':%d:'%(hId)
1415                Limits = calcControls[hfx+'Limits']
1416                x,y,w,yc,yb,yd = Histogram['Data']
1417                xB = np.searchsorted(x,Limits[0])
1418                xF = np.searchsorted(x,Limits[1])
1419                dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
1420                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1421                if len(dMdv):
1422                    dMdv = np.concatenate((dMdv,dMdvh))
1423                else:
1424                    dMdv = dMdvh
1425        return dMdv
1426   
1427    def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
1428        parmdict.update(zip(varylist,values))
1429        G2mv.Dict2Map(parmDict)
1430        Histograms,Phases = HistoPhases
1431        M = np.empty(0)
1432        sumwYo = 0
1433        Nobs = 0
1434        for histogram in Histograms:
1435            if 'PWDR' in histogram[:4]:
1436                Histogram = Histograms[histogram]
1437                hId = Histogram['hId']
1438                hfx = ':%d:'%(hId)
1439                Limits = calcControls[hfx+'Limits']
1440                x,y,w,yc,yb,yd = Histogram['Data']
1441                yc *= 0.0                           #zero full calcd profiles
1442                yb *= 0.0
1443                yd *= 0.0
1444                xB = np.searchsorted(x,Limits[0])
1445                xF = np.searchsorted(x,Limits[1])
1446                Nobs += xF-xB
1447                Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
1448                sumwYo += Histogram['sumwYo']
1449                yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
1450                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1451                yc[xB:xF] *= parmDict[hfx+'Scale']
1452                yc[xB:xF] += yb[xB:xF]
1453                yd[xB:xF] = yc[xB:xF]-y[xB:xF]          #yc-yo then all dydv have no '-' needed
1454                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
1455                M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF])))
1456        Histograms['sumwYo'] = sumwYo
1457        Histograms['Nobs'] = Nobs
1458        Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
1459        if dlg:
1460            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
1461            if not GoOn:
1462                return -M           #abort!!
1463        return M
1464   
1465    ShowBanner()
1466    varyList = []
1467    parmDict = {}
1468    calcControls = {}   
1469    Controls = GetControls(GPXfile)
1470    ShowControls(Controls)           
1471    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
1472    if not Phases:
1473        print ' *** ERROR - you have no histograms to refine! ***'
1474        print ' *** Refine aborted ***'
1475        raise Exception
1476    if not Histograms:
1477        print ' *** ERROR - you have no data to refine with! ***'
1478        print ' *** Refine aborted ***'
1479        raise Exception
1480    phaseVary,phaseDict,pawleyLookup = GetPhaseData(Phases)
1481    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
1482    calcControls.update(controlDict)
1483    histVary,histDict,controlDict = GetHistogramData(Histograms)
1484    calcControls.update(controlDict)
1485    varyList = phaseVary+histVary+hapVary
1486    parmDict.update(phaseDict)
1487    parmDict.update(hapDict)
1488    parmDict.update(histDict)
1489    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
1490    groups,parmlist = G2mv.GroupConstraints(constrDict)
1491    G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
1492    G2mv.Map2Dict(parmDict,varyList)
1493
1494    while True:
1495        begin = time.time()
1496        values =  np.array(Dict2Values(parmDict, varyList))
1497        dlg = wx.ProgressDialog('Residual','Powder profile Rwp =',101.0, 
1498            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
1499        screenSize = wx.ClientDisplayRect()
1500        Size = dlg.GetSize()
1501        dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
1502        Ftol = Controls['min dM/M']
1503        Factor = Controls['shift factor']
1504        try:
1505            if Controls['deriv type'] == 'analytic':
1506                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
1507                    ftol=Ftol,col_deriv=True,factor=Factor,
1508                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1509                ncyc = int(result[2]['nfev']/2)               
1510            else:           #'numeric'
1511                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
1512                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1513                ncyc = int(result[2]['nfev']/len(varyList))
1514        finally:
1515            dlg.Destroy()
1516        runtime = time.time()-begin
1517        chisq = np.sum(result[2]['fvec']**2)
1518        Values2Dict(parmDict, varyList, result[0])
1519        G2mv.Dict2Map(parmDict)
1520        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
1521        GOF = chisq/(Histograms['Nobs']-len(varyList))
1522        print '\n Refinement results:'
1523        print 135*'-'
1524        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
1525        print 'Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1526        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1527        try:
1528            sig = np.sqrt(np.diag(result[1])*GOF)
1529            if np.any(np.isnan(sig)):
1530                print '*** Least squares aborted - some invalid esds possible ***'
1531            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
1532#            for item in table: print item,table[item]               #useful debug - are things shifting?
1533            break                   #refinement succeeded - finish up!
1534        except ValueError:          #result[1] is None on singular matrix
1535            print '**** Refinement failed - singular matrix ****'
1536            Ipvt = result[2]['ipvt']
1537            for i,ipvt in enumerate(Ipvt):
1538                if not np.sum(result[2]['fjac'],axis=1)[i]:
1539                    print 'Removing parameter: ',varyList[ipvt-1]
1540                    del(varyList[ipvt-1])
1541                    break
1542
1543    sigDict = dict(zip(varyList,sig))
1544    SetPhaseData(parmDict,sigDict,Phases)
1545    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
1546    SetHistogramData(parmDict,sigDict,Histograms)
1547    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases)
1548#for testing purposes!!!
1549#    import cPickle
1550#    file = open('structTestdata.dat','wb')
1551#    cPickle.dump(parmDict,file,1)
1552#    cPickle.dump(varyList,file,1)
1553#    for histogram in Histograms:
1554#        if 'PWDR' in histogram[:4]:
1555#            Histogram = Histograms[histogram]
1556#    cPickle.dump(Histogram,file,1)
1557#    cPickle.dump(Phases,file,1)
1558#    cPickle.dump(calcControls,file,1)
1559#    cPickle.dump(pawleyLookup,file,1)
1560#    file.close()
1561
1562def main():
1563    arg = sys.argv
1564    if len(arg) > 1:
1565        GPXfile = arg[1]
1566        if not ospath.exists(GPXfile):
1567            print 'ERROR - ',GPXfile," doesn't exist!"
1568            exit()
1569        GPXpath = ospath.dirname(arg[1])
1570        Refine(GPXfile)
1571    else:
1572        print 'ERROR - missing filename'
1573        exit()
1574    print "Done"
1575         
1576if __name__ == '__main__':
1577    main()
Note: See TracBrowser for help on using the repository browser.