source: trunk/GSASIIstruct.py @ 360

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

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