source: trunk/GSASIIstruct.py @ 369

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

remove wx from GSASIIstruct.py

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