source: trunk/GSASIIstruct.py @ 368

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

fix controls setup & refine enable/disable
add print option to GSASIIstruct.py input routines
needed for new parameter summary feature in GSASII.py
begin view parmater summary stuff

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