source: trunk/GSASIIstruct.py @ 366

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

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