source: trunk/GSASIIstruct.py @ 357

Last change on this file since 357 was 357, checked in by vondreele, 11 years ago

GSASIIpwd.py - add polarization derivative
GSASIIstruct.py - almost complete derivatives for Pawley refinement checked; lacks derivatives for ellipsoidal size parameters

  • 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-08-30 14:38:39 +0000 (Tue, 30 Aug 2011) $
4# $Author: vondreele $
5# $Revision: 357 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 357 2011-08-30 14:38:39Z 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        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
1130        return sig,gam
1131               
1132    hId = Histogram['hId']
1133    hfx = ':%d:'%(hId)
1134    bakType = calcControls[hfx+'bakType']
1135    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1136    yc = np.zeros_like(yb)
1137       
1138    if 'C' in calcControls[hfx+'histType']:   
1139        shl = max(parmDict[hfx+'SH/L'],0.0005)
1140        Ka2 = False
1141        if hfx+'Lam1' in parmDict.keys():
1142            wave = parmDict[hfx+'Lam1']
1143            Ka2 = True
1144            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1145            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1146        else:
1147            wave = parmDict[hfx+'Lam']
1148    else:
1149        print 'TOF Undefined at present'
1150        raise ValueError
1151    for phase in Histogram['Reflection Lists']:
1152        refList = Histogram['Reflection Lists'][phase]
1153        Phase = Phases[phase]
1154        pId = Phase['pId']
1155        pfx = '%d::'%(pId)
1156        phfx = '%d:%d:'%(pId,hId)
1157        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1158        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1159        sizeEllipse = []
1160        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1161            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1162        for refl in refList:
1163            if 'C' in calcControls[hfx+'histType']:
1164                h,k,l = refl[:3]
1165                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1166                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
1167                Icorr = GetIntensityCorr(refl,phfx,hfx,calcControls,parmDict)
1168                if 'Pawley' in Phase['General']['Type']:
1169                    try:
1170                        refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
1171                    except KeyError:
1172                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1173                        continue
1174                else:
1175                    raise ValueError       #wants strctrfacr calc here
1176                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1177                iBeg = np.searchsorted(x,refl[5]-fmin)
1178                iFin = np.searchsorted(x,refl[5]+fmax)
1179                if not iBeg+iFin:       #peak below low limit - skip peak
1180                    continue
1181                elif not iBeg-iFin:     #peak above high limit - done
1182                    return yc,yb
1183                yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1184                if Ka2:
1185                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1186                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
1187                    iBeg = np.searchsorted(x,pos2-fmin)
1188                    iFin = np.searchsorted(x,pos2+fmax)
1189                    yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
1190            else:
1191                raise ValueError
1192    return yc,yb   
1193           
1194def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1195   
1196    def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1197        gamDict = {}
1198        costh = cosd(refl[5]/2.)
1199        tanth = tand(refl[5]/2.)
1200        if calcControls[phfx+'SizeType'] == 'isotropic':
1201            gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh)
1202            gamDict[phfx+'Size:0'] = -gam/parmDict[phfx+'Size:0']
1203        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1204            H = np.array(refl[:3])
1205            P = np.array(calcControls[phfx+'SizeAxis'])
1206            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1207            Si = parmDict[phfx+'Size:0']
1208            Sa = parmDict[phfx+'Size:1']
1209            gami = (1.8*wave/np.pi)/(Si*Sa)
1210            sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2)
1211            gam = gami*sqtrm/costh           
1212            gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)-gam/Si
1213            gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)-gam/Sa         
1214        else:           #ellipsoidal crystallites - do numerically?
1215            H = np.array(refl[:3])
1216            gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1217                       
1218        if calcControls[phfx+'MustrainType'] == 'isotropic':
1219            gamDict[phfx+'Mustrain:0'] =  0.018*tanth/np.pi           
1220        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1221            H = np.array(refl[:3])
1222            P = np.array(calcControls[phfx+'MustrainAxis'])
1223            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1224            Si = parmDict[phfx+'Mustrain:0']
1225            Sa = parmDict[phfx+'Mustrain:1']
1226            gami = 0.018*Si*Sa*tanth/np.pi
1227            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1228            gam = gami/sqtrm
1229            gamDict[phfx+'Mustrain:0'] = gam/Si-gami*Si*cosP**2/sqtrm**3
1230            gamDict[phfx+'Mustrain:1'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
1231        else:       #generalized - P.W. Stephens model
1232            pwrs = calcControls[phfx+'MuPwrs']
1233            const = 0.018*refl[4]**2*tanth
1234            for i,pwr in enumerate(pwrs):
1235                gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1236        return gamDict
1237       
1238    def GetIntensityDerv(refl,phfx,hfx,calcControls,parmDict):
1239        Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1240        pola,dpdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1241        dIdpola = Icorr*dpdPola
1242        Icorr *= pola
1243        dIdsh = Icorr/parmDict[hfx+'Scale']
1244        dIdsp = Icorr/parmDict[phfx+'Scale']
1245       
1246        return Icorr,dIdsh,dIdsp,dIdpola
1247       
1248    def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1249        dpr = 180./np.pi
1250        h,k,l = refl[:3]
1251        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1252        dst = np.sqrt(dstsq)
1253        pos = refl[5]
1254        const = dpr/np.sqrt(1.0-wave*dst/4.0)
1255        dpdw = const*dst
1256        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1257        dpdA *= const*wave/(2.0*dst)
1258        dpdZ = 1.0
1259        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1260        if 'Bragg' in calcControls[hfx+'instType']:
1261            dpdSh = -4.*const*cosd(pos/2.0)
1262            dpdTr = -1.e-7*const*sind(pos)
1263            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1264        else:               #Debye-Scherrer - simple but maybe not right
1265            dpdXd = -const*cosd(pos)
1266            dpdYd = -const*sind(pos)
1267            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1268           
1269    def cellVaryDerv(pfx,SGData,dpdA): 
1270        if SGData['SGLaue'] in ['-1',]:
1271            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1272                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1273        elif SGData['SGLaue'] in ['2/m',]:
1274            if SGData['SGUniq'] == 'a':
1275                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1276            elif SGData['SGUniq'] == 'b':
1277                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1278            else:
1279                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1280        elif SGData['SGLaue'] in ['mmm',]:
1281            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1282        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1283            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
1284        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1285            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
1286        elif SGData['SGLaue'] in ['3R', '3mR']:
1287            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1288        elif SGData['SGLaue'] in ['m3m','m3']:
1289            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
1290   
1291    lenX = len(x)               
1292    hId = Histogram['hId']
1293    hfx = ':%d:'%(hId)
1294    bakType = calcControls[hfx+'bakType']
1295    dMdv = np.zeros(shape=(len(varylist),len(x)))
1296    if hfx+'Back:0' in varylist:
1297        dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1298        bBpos =varylist.index(hfx+'Back:0')
1299        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1300       
1301    if 'C' in calcControls[hfx+'histType']:   
1302        dx = x[1]-x[0]
1303        shl = max(parmDict[hfx+'SH/L'],0.002)
1304        Ka2 = False
1305        if hfx+'Lam1' in parmDict.keys():
1306            wave = parmDict[hfx+'Lam1']
1307            Ka2 = True
1308            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1309            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1310        else:
1311            wave = parmDict[hfx+'Lam']
1312    else:
1313        print 'TOF Undefined at present'
1314        raise ValueError
1315    for phase in Histogram['Reflection Lists']:
1316        refList = Histogram['Reflection Lists'][phase]
1317        Phase = Phases[phase]
1318        SGData = Phase['General']['SGData']
1319        pId = Phase['pId']
1320        pfx = '%d::'%(pId)
1321        phfx = '%d:%d:'%(pId,hId)
1322        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1323        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1324        sizeEllipse = []
1325        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1326            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1327        for iref,refl in enumerate(refList):
1328            if 'C' in calcControls[hfx+'histType']:
1329                Icorr,dIdsh,dIdsp,dIdpola = GetIntensityDerv(refl,phfx,hfx,calcControls,parmDict)
1330                hkl = refl[:3]
1331                pos = refl[5]
1332                tanth = tand(pos/2.0)
1333                costh = cosd(pos/2.0)
1334                dsdU = tanth**2
1335                dsdV = tanth
1336                dsdW = 1.0
1337                dgdX = 1.0/costh
1338                dgdY = tanth
1339                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1340                iBeg = np.searchsorted(x,refl[5]-fmin)
1341                iFin = np.searchsorted(x,refl[5]+fmax)
1342                dMdpk = np.zeros(shape=(6,len(x)))
1343                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
1344                for i in range(1,5):
1345                    dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i]
1346                dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*dMdipk[0]
1347                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1348                if Ka2:
1349                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1350                    kdelt = int((pos2-refl[5])/dx)               
1351                    iBeg = min(lenX,iBeg+kdelt)
1352                    iFin = min(lenX,iFin+kdelt)
1353                    if iBeg-iFin:
1354                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,xdata[iBeg:iFin])
1355                        for i in range(1,5):
1356                            dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]
1357                        dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
1358                        dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
1359                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*Icorr*refl[8]}
1360                try:
1361                    idx = varylist.index(pfx+'PWLref:'+str(iref))
1362                    dMdv[idx] = dervDict['int']
1363                except ValueError:
1364                    pass
1365                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1366                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1367                    hfx+'U':[dsdU,'sig'],hfx+'V':[dsdV,'sig'],hfx+'W':[dsdW,'sig'],
1368                    hfx+'X':[dgdX,'gam'],hfx+'Y':[dgdY,'gam'],hfx+'SH/L':[1.0,'shl'],
1369                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1370                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1371                    hfx+'DisplaceY':[dpdY,'pos'],}
1372                for name in names:
1373                    if name in varylist:
1374                        item = names[name]
1375                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
1376                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1377                for name,dpdA in cellDervNames:
1378                    if name in varylist:
1379                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
1380                gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse)
1381                for name in gamDict:
1382                    if name in varylist:
1383                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']                       
1384            else:
1385                raise ValueError
1386           
1387    return dMdv   
1388                   
1389def Refine(GPXfile):
1390   
1391    def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
1392        parmdict.update(zip(varylist,values))
1393        G2mv.Dict2Map(parmDict)
1394        Histograms,Phases = HistoPhases
1395        dMdv = np.empty(0)
1396        for histogram in Histograms:
1397            if 'PWDR' in histogram[:4]:
1398                Histogram = Histograms[histogram]
1399                hId = Histogram['hId']
1400                hfx = ':%d:'%(hId)
1401                Limits = calcControls[hfx+'Limits']
1402                x,y,w,yc,yb,yd = Histogram['Data']
1403                xB = np.searchsorted(x,Limits[0])
1404                xF = np.searchsorted(x,Limits[1])
1405                dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
1406                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1407                if len(dMdv):
1408                    dMdv = np.concatenate((dMdv,dMdvh))
1409                else:
1410                    dMdv = dMdvh
1411        return dMdv
1412   
1413    def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
1414        parmdict.update(zip(varylist,values))
1415        G2mv.Dict2Map(parmDict)
1416        Histograms,Phases = HistoPhases
1417        M = np.empty(0)
1418        sumwYo = 0
1419        Nobs = 0
1420        for histogram in Histograms:
1421            if 'PWDR' in histogram[:4]:
1422                Histogram = Histograms[histogram]
1423                hId = Histogram['hId']
1424                hfx = ':%d:'%(hId)
1425                Limits = calcControls[hfx+'Limits']
1426                x,y,w,yc,yb,yd = Histogram['Data']
1427                yc *= 0.0                           #zero full calcd profiles
1428                yb *= 0.0
1429                yd *= 0.0
1430                xB = np.searchsorted(x,Limits[0])
1431                xF = np.searchsorted(x,Limits[1])
1432                Nobs += xF-xB
1433                Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
1434                sumwYo += Histogram['sumwYo']
1435                yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
1436                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1437                yc[xB:xF] *= parmDict[hfx+'Scale']
1438                yc[xB:xF] += yb[xB:xF]
1439                yd[xB:xF] = yc[xB:xF]-y[xB:xF]          #yc-yo then all dydv have no '-' needed
1440                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
1441                M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF])))
1442        Histograms['sumwYo'] = sumwYo
1443        Histograms['Nobs'] = Nobs
1444        Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
1445        if dlg:
1446            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
1447            if not GoOn:
1448                return -M           #abort!!
1449        return M
1450   
1451    ShowBanner()
1452    varyList = []
1453    parmDict = {}
1454    calcControls = {}   
1455    Controls = GetControls(GPXfile)
1456    ShowControls(Controls)           
1457    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
1458    if not Phases:
1459        print ' *** ERROR - you have no histograms to refine! ***'
1460        print ' *** Refine aborted ***'
1461        raise Exception
1462    if not Histograms:
1463        print ' *** ERROR - you have no data to refine with! ***'
1464        print ' *** Refine aborted ***'
1465        raise Exception
1466    phaseVary,phaseDict,pawleyLookup = GetPhaseData(Phases)
1467    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
1468    calcControls.update(controlDict)
1469    histVary,histDict,controlDict = GetHistogramData(Histograms)
1470    calcControls.update(controlDict)
1471    varyList = phaseVary+histVary+hapVary
1472    parmDict.update(phaseDict)
1473    parmDict.update(hapDict)
1474    parmDict.update(histDict)
1475    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
1476    groups,parmlist = G2mv.GroupConstraints(constrDict)
1477    G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
1478    G2mv.Map2Dict(parmDict,varyList)
1479
1480    while True:
1481        begin = time.time()
1482        values =  np.array(Dict2Values(parmDict, varyList))
1483        dlg = wx.ProgressDialog('Residual','Powder profile Rwp =',101.0, 
1484            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
1485        screenSize = wx.ClientDisplayRect()
1486        Size = dlg.GetSize()
1487        dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
1488        Ftol = Controls['min dM/M']
1489        try:
1490            if Controls['deriv type'] == 'analytic':
1491                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,ftol=Ftol,col_deriv=True,
1492                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1493                ncyc = int(result[2]['nfev']/2)               
1494            else:           #'numeric'
1495                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,
1496                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1497                ncyc = int(result[2]['nfev']/len(varyList))
1498        finally:
1499            dlg.Destroy()
1500        runtime = time.time()-begin
1501        chisq = np.sum(result[2]['fvec']**2)
1502        Values2Dict(parmDict, varyList, result[0])
1503        G2mv.Dict2Map(parmDict)
1504        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
1505        GOF = chisq/(Histograms['Nobs']-len(varyList))
1506        print '\n Refinement results:'
1507        print 135*'-'
1508        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
1509        print 'Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1510        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1511        try:
1512            sig = np.sqrt(np.diag(result[1])*GOF)
1513            if np.any(np.isnan(sig)):
1514                print '*** Least squares aborted - some invalid esds possible ***'
1515            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
1516#            for item in table: print item,table[item]               #useful debug - are things shifting?
1517            break                   #refinement succeeded - finish up!
1518        except ValueError:          #result[1] is None on singular matrix
1519            print '**** Refinement failed - singular matrix ****'
1520            Ipvt = result[2]['ipvt']
1521            for i,ipvt in enumerate(Ipvt):
1522                if not np.sum(result[2]['fjac'],axis=1)[i]:
1523                    print 'Removing parameter: ',varyList[ipvt-1]
1524                    del(varyList[ipvt-1])
1525                    break
1526
1527    sigDict = dict(zip(varyList,sig))
1528    SetPhaseData(parmDict,sigDict,Phases)
1529    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
1530    SetHistogramData(parmDict,sigDict,Histograms)
1531    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases)
1532#for testing purposes!!!
1533#    import cPickle
1534#    file = open('structTestdata.dat','wb')
1535#    cPickle.dump(parmDict,file,1)
1536#    cPickle.dump(varyList,file,1)
1537#    for histogram in Histograms:
1538#        if 'PWDR' in histogram[:4]:
1539#            Histogram = Histograms[histogram]
1540#    cPickle.dump(Histogram,file,1)
1541#    cPickle.dump(Phases,file,1)
1542#    cPickle.dump(calcControls,file,1)
1543#    cPickle.dump(pawleyLookup,file,1)
1544#    file.close()
1545
1546def main():
1547    arg = sys.argv
1548    if len(arg) > 1:
1549        GPXfile = arg[1]
1550        if not ospath.exists(GPXfile):
1551            print 'ERROR - ',GPXfile," doesn't exist!"
1552            exit()
1553        GPXpath = ospath.dirname(arg[1])
1554        Refine(GPXfile)
1555    else:
1556        print 'ERROR - missing filename'
1557        exit()
1558    print "Done"
1559         
1560if __name__ == '__main__':
1561    main()
Note: See TracBrowser for help on using the repository browser.